Electronic and vibrational spectroscopy of miscible MgO–ZnO ternary alloys

Zinc oxide (ZnO) is a wide, direct bandgap $Eg∼$3.4 eV II–VI semiconductor that has attracted increasing interest for both fundamental research and diverse optoelectronic applications.^1–3 For example, ZnO is a promising material for light-emitting and laser diodes operating in blue and ultraviolet (UV) regimes.^2–4 However, low electroluminescence often diminishes the efficiency of ZnO-based devices.^5,6 To realize and enhance the efficiency of ZnO-based devices, their electronic structure, through band-engineering, can be systematically tuned by doping. ZnO can be alloyed with MgO to form a ternary $MgxZn1−xO$ compound, thus, enabling bandgap engineering and luminescence in the UV regime.⁷ Often, the combination of these group-II oxides in alloys leads to crystal structure mismatch: undoped ZnO prefers the hexagonal wurtzite (B4) structure or the fourfold coordinated zinc-blende structure (B3), while MgO favors the cubic rock salt (B1) structure at ambient conditions.^8,9 Experiments have shown that $MgxZn1−xO$ exhibits the B4 structure for high ZnO concentration,¹⁰ while preferring the B1 structure at high MgO concentration.¹¹ At intermediate concentrations, it exhibits phase separation¹² due to its compositional gradient, which often leads to a thermodynamically unstable crystal structure. In general, the isovalent and isostructural II–VI alloys are thermodynamically unstable because the mixing enthalpy in either the B1, B3, or B4 structure is always positive.^12,14

The difficulty in synthesizing stable phases of MgO–ZnO alloys has hindered the adequate exploitation of their promising properties for technological applications. This prompted earlier calculations on the stability of the $MgxZn1−xO$ alloy system.^12,13 For example, the work of Sanati et al.¹² on various alloys of MgO–ZnO found energetic stability, under certain conditions, in the sixfold-coordinated structure for Zn concentrations below 67%. Instead of the often-used traditional doping approach of alloying ZnO with MgO to obtain the phase immiscible Mg_xZn_1−xO alloy, we can apply high-throughput crystal structure prediction techniques to search for energetically and dynamically stable stoichiometric MgO–ZnO alloy systems. A high-throughput experiment to achieve this is, however, daunting because of the ample parameter space that needs to be explored. Computationally, such an approach is now routinely used to predict new crystal structures,^15,16 including, recently, the prediction of energetically stable high pressure phase of $MgxZn1−xO$⁠.¹⁷ Herein, we report the electronic structure and thermodynamic stability of the stable phases of MgZnO$2$ crystal structure at ambient conditions obtained from the recently developed evolutionary algorithm for predicting crystal structure.¹⁶ We provide detailed physical properties and characterize electronic, mechanical, phononic, and optical spectroscopy. We extensively discuss the results in relation to experimental data of the parent systems and provide details of the crystal structure information (see supplementary material). We hope that they will motivate future experimental investigations of the electronic structure, particularly using optical measurements and photoemission, and vibrational spectroscopy by Raman measurements.

The ground state crystal structures of MgZnO$2$ are predicted using the ab initio evolutionary algorithm CALYPSO.¹⁶ We performed variable compositional simulations for $MgnZnmOm+n$ systems (n, m $≤30$⁠) while maintaining the stoichiometric composition of MgZnO$2$⁠. The structural optimizations and electronic structure calculations for the selected structures are performed using density functional theory (DFT) as implemented in the Vienna Ab Initio Simulation Package (VASP) code.¹⁸ The electronic wavefunctions of the systems are represented by a plane wave basis with an energy cutoff of 550 eV. The exchange-correlation interactions were described with the HSE06 hybrid functional¹⁹ with a $10×10×10$ $Γ$—centered $k$—point grid, which ensures that the energy (charge) is converged to within 10$−4$ (10$−8$⁠) eV. We also obtained the electronic structure using the strongly constrained and appropriately normed (SCAN) semilocal density functional²⁰ (see the supplementary material). The phononic and vibrational spectroscopy is obtained using density functional perturbation theory (DFPT)²¹ as implemented in VASP in conjunction with the Phonopy package.²² 

In this work, we present two stable ground state crystal structures of MgZnO$2$⁠. The two of the predicted structures belong to the tetragonal crystal symmetry $I41/amd$ space group No. 141 and P$4¯$m2 space group No. 115, respectively (see, Fig. 1). We designate the former and latter structures as A and B, respectively. The predicted lattice constants, along with other salient parameters, are presented in Table I, and the Wyckoff positions in the crystallographic information file format are included in the supplementary material. The phase stability of the MgO–ZnO systems is studied by calculating the formation energy $Eform$ per atom; $Eform$ is a measure of the energetic stability of a material and for our system, it is obtained as $Eform=Etot−∑iciμi$⁠, where $i$ denotes different types of atoms present in the unit cell of the material with concentration $ci$ and $μi$ is the chemical potential approximated with the standard state (bulk) energy $Ei$ of the corresponding atom, $i$⁠. The two structures exhibit negative formation energies with an average $Eform∼−2.11$eV/atom; this implies an exothermic process confirming the phase and energetic stability of the crystal. To further confirm the energetic stability for the predicted structures, we calculate the enthalpy profile as $ΔH=ΔG+TΔS$⁠, where $H$⁠, $G$⁠, $S$⁠, and $T$ is enthalpy, free energy, entropy, and temperature, respectively, for DFT, $T=0$K as a function of generations (see Fig. S1 in the supplementary material). All the predicted structures show a negative formation enthalpy $ΔH>−5.0$ eV, which supports the energetic stability of the structures. The predicted $ΔH$ is consistent with previous theoretical results for the Mg$1−x$Zn$x$O system;^12,13 though these previously reported $ΔH$ are negative, they are smaller than our results due to different stoichiometric compositions.

While the energetic stability is often a confirmation that a crystal can be experimentally synthesized, in most cases, basing the overall stability on the energy alone may be misleading. In the absence of any phase transition, a crystal may have negative formation energy yet thermodynamically unstable. To ascertain the dynamical stability of the predicted structures, we performed phonon and vibrational spectroscopy calculations (Fig. 2). We also show in Fig. 2 the Raman and Infrared spectroscopy. Both structures show positive definite phonon frequency along the various high symmetry point in the first Brillouin zone, confirming the dynamical stability. From the phonon atomic projected density of states, the acoustic branch is dominated by Zn atoms, while the optical branch is composed of a strong hybridization between the O and Mg atoms. This trend is consistent with the domination of the frequency scale of the acoustic (optical) phonon modes by atoms with larger (smaller) masses similar to a diatomic linear chain model.²³ 

The predicted crystal structures belong to the tetragonal crystal symmetry, albeit with a different point group symmetry leading to a different irreducible representation of the vibrational modes. Crystal structure A belongs to the point group D$4h$ (4/mmm) with 16 atoms in the primitive cell. The irreducible representation of the corresponding 48 normal vibrational modes at the Brillouin zone center (⁠$Γ$-point) is $Γ=2A1g(R)+3A1u+A2g+6A2u(IR)+2B1g(R)+3B1u+B2g(R)+6B2u+9Eu(IR)+3Eg(R)$⁠, where IR and R are the infrared and Raman active modes. $A1u$⁠, $A2g$⁠, $B1u$⁠, and $B2u$ are the silent modes. The R-active $Eg$⁠, $B1g$⁠, and the $A1g$ modes are mainly due to the symmetric stretching vibration, the symmetric bending vibration, and the antisymmetric bending vibration of the O–Mg/Zn–O, respectively. The phonon spectra, the phonon total, and the atomic projected density of states along with the Raman and Infrared spectra are presented in Fig. 2(a). (Further analysis of the irreducible representation of the phonon modes are shown in Table S1 of the supplementary material.) We hope the presented phonon, vibrational, Raman, and Infrared spectroscopy will guide experimental and other analyses. For the crystal structure B with point group symmetry $D2d(−42m)$⁠, the primitive cell contains four atoms; as such, there will be a total of twelve normal vibrational modes at the $Γ$ point [Fig. 2(b)]. The irreducible representation of the vibrational modes at the $Γ$ point could be expressed as $Γ=A1(R)+3B2(IR+R)+4E(IR+R)$⁠. One low-energy $B2$ and one double $E$ modes are the usual acoustic modes characterized by the transverse acoustic (TA), longitudinal acoustic (LA), and the out-of-plane transverse acoustic or flexural mode (ZA); the remaining belong to the optical mode. Using group theory analysis²⁴ in conjunction with the Raman and Infrared spectra [Fig. 2(b)], we assign the feature at $∼178.70$ cm$−1$ to the degenerate IR- and R-active $E$ mode and the structure at $∼326.28$ cm$−1$ to the IR- and R-active $B2$ mode. We predict a degenerate IR- and R-active $E$ mode at $∼413.94$ cm$−1$⁠, followed by IR- and R-active $E$ mode at $∼480.60$ cm$−1$⁠. We predict an R-active $A1$ mode though with small activity around 531.00 cm$−1$ before a degenerate IR- and R-active $E$ mode at $∼548.81$ cm$−1$⁠. (Further details of the irreducible representation are presented in Table S2 in the supplementary material.) We note that the predicted features in both the vibrational and Raman spectra depict more of ZnO character than MgO. For example, room temperature Raman spectroscopy experiments reported the R-active modes $E∼583$ cm$−1$ and $A1∼574$ cm$−1$⁠.²⁵ 

The Debye temperature $ΘD$ is an important thermal property of a material. It quantifies the temperature below which modes start to freeze out and above which, the modes gain enough energy to be excited.²⁶ The Debye temperature can be obtained as $ΘD3=13∑i1ΘΓ,i3$⁠, where $ΘΓ,i=ℏΩΓ,i/κB$ and $ΘΓ,i$ corresponds to the frequency at the zone boundary of the $i$th acoustic mode.²⁷ The predicted $ΘD$ is presented in Table I, which seem to be closer to that of ZnO with $ΘD∼440$ K than to MgO with $ΘD∼948$ K as reported by experiments.^9,25

The elastic properties not only provide the information on the mechanical stability of material but also serve as a guide on the potential device applications of a given material. In this regard, we have calculated the elastic properties of the studied materials using the finite difference method. We used a plane wave energy cutoff of 700 eV to ensure the adequate convergence of the stress tensor. The elastic tensor matrix is determined by performing six finite distortions of the lattice and obtaining the elastic constants from the strain–stress relationship.^28,29 The predicted elastic tensor matrices are presented in Table II. Using the calculated elastic tensor matrix, we predict the elastic constants presented in Table I. Aside from the isotropic elastic behavior [the calculated elastic anisotropy ratio $A=2C44/(C11−C12)∼1.0$], the two studied structures show markedly different elastic behavior. Structure A (B) could withstand linear compressibility of $∼$1.81 (1.65 TPa$−1$⁠). The predicted Bulk modulus of structure A is $∼41$ GPa larger than that of structure B. The same trend is observed for Young’s modulus E, shear modulus G, P-wave modulus M, and the Poisson ratio $μ$⁠, which is larger by $E∼107.37$ GPa, $G∼43.65$ GPa, $M∼98.89$ GPa, and $μ∼0.07$⁠. Though structure A exhibits large elastic constants, structure B is, however, more ductile from the bulk/shear ratio value (Table I).

The description of the electronic properties of the two forms of MgZnO$2$ studied herein as obtained using the HSE06 hybrid functional is provided by the band structure, total, and atomic projected density of states in Fig. 3. We predict an indirect energy bandgap $Eg∼4.62$ eV along the $Γ−N$ point of the high symmetry point of the Brillouin zone for structure A. We note that the direct energy bandgap is 4.52 eV, which is only 0.10 eV smaller than the largest bandgap; as such, the indirect and direct energy bandgaps could be said to be almost degenerate. Structure B is predicted to be a direct bandgap semiconductor with $Eg∼3.70$ eV at the $Γ$-point. The predicated bandgaps are listed in Table I. We note that the energy bandgap predicted with the semilocal SCAN (see Fig. S2 in the supplementary material) and the Perdew–Burke–Ernzerhof generalized gradient approximation (PBE-GGA)³⁰ is significantly smaller in both materials. For example, for structure A, SCAN and PBE-GGA functional led to $Eg∼3.32$ and 2.73 eV, respectively. We obtained $Eg∼2.43$ and 1.98 eV using SCAN and PBE-GGA functional, respectively, for structure B.

The two studied structures share in common some of the features of ZnO. For example, the predicted energy bandgap is in the same range as the $∼3.4$ eV reported for bulk ZnO by several experiments and computations.^9,25,31–33 Aside from the difference in the energy bandgap of the two structures, the character of the states in the proximity of the Fermi energy are similar for the two structures. From the projected density of states, the lower-lying group in the valence states consists of highly hybridized O (from s- and p-states), while the upper group is mostly dominated by Zn (d-states), O (p-states), and tiny contribution from Mg (s- and p-states). The uppermost conduction states are dominated by a strong hybridization between Zn (s-states) and O (s- and p-states). This is evident from the free electron-like (s-like) states at the conduction band minimum (CBM) located at the $Γ$-point. The total width of the upper valence bands within 10.0 eV from the valence band maximum (VBM) is 6.30 (5.28 eV) for structure A (B). This value is in good agreement with $∼6.50$ eV for MgO and approximately 1–2 eV smaller than $∼7.2$ eV for ZnO reported by experiments.^31,34 The total width of the entire valence states is 19.52 (18.25 eV) for structure A (B). These values are within the experimental uncertainty of the total valence width of ZnO.²⁵ We note that the valence bandwidth obtained with the scan functional decreased by $∼$5.2%. With the PBE-GGA functional, the decrease is more significant $∼14.07$%.

The effective mass is a measure of the curvature of the calculated band structure, and it provides insight into the transport properties of materials. To this end, we calculate the carrier band effective mass of the two studied structures. The band effective mass $mb$ is obtained from the band structure by fitting a parabola $Ek=ℏ22m0k→TAk→$ to the states around the band extremum (CBM and VBM) of Fig. 3, where $k=(∥,⊥)$ is the parallel and perpendicular direction of the $k$-point measured from the band extremum and $m0$ is the free electron mass. The predicted electron and hole effective masses along the $∥$ (⁠$⊥$⁠) direction are listed in Table I. The effective mass decreased by $∼6.50$ and 15.20%, respectively, with the SCAN and PBE-GGA functionals. The large hole effective mass, especially for structure A, is due to the rather flatbands around the top of the valence states. We note that the predicted $me∥/⊥$ straddle that of ZnO $∼0.275$ and MgO $∼0.370$ m$0$⁠.^9,25,35,36 While the predicted $mh∥$ is in good agreement with $∼0.60$ m$0$ reported for ZnO,²⁵ the calculated $mh⊥$ is in better agreement with the 1–2 m$0$ reported for MgO.³⁶ 

To gain further insights into the physical properties of the predicted structures, we calculate the electrodynamical response as a function of photon energy (⁠$ℏω$⁠). The optical spectroscopy provides a detailed and direct information on charge dynamics and “bulk” electronic structure of materials since it probes both free carriers and interband excitation. The optical response is determined using the spectra from our electronic structure calculations obtained with HSE06 hybrid functional without any shift or scissor operator. We used a small broadening parameter $∼10−4$ eV as such, our optical spectra show a lot of features. Using random phase approximation, we compute the real part of the complex dynamical dielectric function

where $Ξ=⟨uck+eαq|uvk⟩⟨uck+eβq|uvk⟩∗$⁠, $ω$ is the photon frequency, $α$ (⁠$β$⁠) is band indices, $v$ and $c$ denote valence and conduction band states, and $uck$ is the cell periodicity at the $k$-point over volume $Ω$⁠. The dispersive (real) part of the dynamical dielectric function $ϵ(1)(ω)$ is then obtained using the appropriated Kramers–Kronig transformation

where $P$ is the principal value.²¹ The calculated optical spectra are presented in Fig. 4. The onset of absorption that is determined by the absorption edge depicts the characteristic optical bandgap. This corresponds to the lowest sharp structure around 4.53 (3.12) eV for structure A (B) due to the direct optical excitation.

The predicted optical bandgap is within the numerical uncertainty of the ones obtained from our band structure calculations (Table I). The optical spectra for both materials share almost the same characteristic features. In the $ϵ2(ω)$ spectra [Fig. 4(a)], we found a sharp structure around 6.50 (5.06), followed by a set of close sharp features and shoulders that transcend in both structures up to 17.5 eV, before the spectra systematically decay toward zero photon energy.

In Fig. 4(b), we present the calculated dispersive parts of the dynamical dielectric function $ϵ(1)(ω)$⁠. The main features in both structures are sharp structures around 4.95 (3.50) eV for structure A (B), series of other features between $∼5.20$ eV and 15.10 eV, characterized by a steep decrease that started around 7.80 eV. Then, $ϵ(1)(ω)$ for both structures becomes negative with a minimum at $∼15.10$ eV. The optical spectra obtained using the electronic structure spectra from the SCAN semilocal functional though exhibits the same trend, the critical energies are significantly different (see Fig. S3 of the supplementary material). We attribute this to the nonlocal effects that are inherent in the HSE06 hybrid functional as evident from the sharper structures in the optical spectra. We further determined the high photon energy $ϵ∞$ and the static photon energy $ϵ0$ limits of the dielectric constants. $ϵ∞$ is obtained as the zero-photon energy limit of the $ϵ(1)(ω)$ spectra and the obtained values are listed in Table I. To obtain $ϵ0$⁠, we note that $ϵ0=ϵ∞+ϵion$⁠, where $ϵion$ is the static photon energy (⁠$ℏω=0$⁠) contribution of the ionic lattice dynamics to the dielectric constant. We obtained from our phonon calculations $ϵion(∥/⊥)∼5.04/4.91$ (3.97/2.57) for structure A (B). Using the obtained $ϵion$⁠, we predict the static photon energy of the dielectric constant summarized in Table I. The predicted optical behavior of structure A is generally isotropic while that of structure B could be said to be slightly anisotropic, mainly from the dynamics of $ϵion(∥/⊥)$⁠. We note that the $ϵ∞(∥/⊥)$ obtained using the eigenstates from the SCAN functional as the input in our optical calculations is 3.43/3.43 (3.34/3.31) for structure A (B), which seem to be closer to the values reported for ZnO.²⁵ Overall, compared to the parent materials, the predicted limits of the dielectric constants straddle those of ZnO and MgO. For example, experiments report ZnO has $ϵ∞(∥/⊥)∼3.75/3.70$ and $ϵ0(∥/⊥)∼8.75/7.80$⁠, while MgO exhibits $ϵ∞(∥/⊥)∼2.94$ and $ϵ0(∥/⊥)∼9.83$⁠.²⁵ 

We report detailed ab initio results of the structural, phononic, electronic, transport, and optical properties of newly predicted crystal structures of the MgZnO$2$ alloy system. Through the calculations of formation energies and the phonon and vibrational spectra, including the Raman and Infrared spectroscopy, we establish the dynamical stability of the predicted structures. We demonstrate highly tunable electronic and optical properties of two of the predicted structures with optoelectronic features that make them promising for application as a photonic material. We provide detailed crystal structure information to guide further experimental studies.

See the supplementary material for the predicted crystal lattice parameters, the irreducible representation of the phonon active modes at the $Γ$-point, and the predicted electronic and optical properties obtained with the SCAN functional.

Supercomputer and computational resources were provided by the Lehigh University (LU) HPC Center. K.A. acknowledges support of the LU Physics REU program (No. PHY-1852010) and C.E.E. acknowledges the LU start-up and summer research fellowship.