First-principles surface energies for monoclinic Ga₂O₃ and Al₂O₃ and consequences for cracking of (Al_xGa_1−x)₂O₃

Due to its large bandgap (4.76 eV–5.04 eV),^1–4 monoclinic β-Ga₂O₃ provides a large critical electric field strength in field effect transistors (FETs).^5,6 Despite being a wide bandgap material, Ga₂O₃ can easily be n-type doped.^5,7 β-Ga₂O₃ thus constitutes a promising platform for applications in high-power electronic devices.

Alloying can be used to tune device performance through bandgap engineering. Alloys with Al have increased bandgaps and enable carrier confinement in heterojunction FETs.^8,9 The smaller Al atomic size leads to smaller lattice parameters in (Al_xGa_1−x)₂O₃. For (Al_xGa_1−x)₂O₃ epilayers pseudomorphically grown on β-Ga₂O₃ substrates with the same crystallographic orientation as the substrate, tensile stress in (Al_xGa_1−x)₂O₃ is expected due to the lattice mismatch.¹⁰ For a thin layer, the stress can be accommodated by elastic deformation. As the elastic energy in the growing film increases, it can lead to the formation of misfit dislocations, surface roughness,¹¹ V-shaped defects,¹¹ or cracks,¹² leading to a critical thickness where the strain can no longer be accommodated elastically. For the growth of β-(Al_xGa_1−x)₂O₃ films on Ga₂O₃ substrates, cracks have been found to be the major limitation.¹³ Understanding the mechanism of crack formation is essential for growing the high-quality films required for devices.

In this paper, we study the critical thickness for (Al_xGa_1−x)₂O₃ films grown on β-Ga₂O₃ substrates along [100], [010], and [001] directions; among these, the [010] orientation is most commonly used in growth. Our approach is based on the energy balance between the elastic energy and the brittle fracture toughness, where the latter is derived from surface energies. All quantities are calculated with first-principles density functional theory (DFT). DFT calculations for surface energies have been previously reported,^14,15 but they did not cover all the orientations discussed here and were not at the same level of accuracy as the results presented here.

The calculations of surface energies performed in the context of crack formation also prompted us to perform an in-depth analysis of surface stability. We explore the mechanisms that determine the relative stability of different surfaces; we analyze bond strengths, dangling-bond densities, and atomic relaxation, providing a comprehensive picture of how these factors impact the energetics. These insights are important for surface science in general, particularly in the context of understanding growth.

We find that the (100)B surface¹⁴ is most stable and acts as the crack plane for (Al_xGa_1−x)₂O₃ film growth along the [010] and [001] directions. The remarkably small energies of (100)B crack planes greatly limit the film thickness and hinder the incorporation of larger concentrations of Al. Growth along the [100] direction gives rise to a considerably larger critical thickness, owing to the higher energy of the favored crack plane, which is (001)B.

Equation (1) represents an energy balance between the elastic energy stored in the strained film per unit area (⁠$h2∑i,jϵiCijϵj$⁠) and the fracture toughness, Γ. The fracture toughness consists of the surface energy when the crack is created as well as the energy of plastic deformation at the crack tip.¹² In this study, we assume brittle fracture (Γ = Γ_b), meaning that only the surface energy is taken into account in the fracture toughness. Therefore, the calculated h_c constitutes a lower bound on the true critical thickness.

First-principles calculations were performed using DFT and projector augmented wave (PAW) potentials¹⁷ as implemented in the Vienna ab initio Simulation Package (VASP).^18,19 A 500 eV kinetic energy cutoff was chosen for the plane wave expansion. The PAW potentials correspond to the valence-electron configurations 3d¹⁰4s²4p¹ for Ga, 3s²3p¹ for Al, and 2s²2p⁴ for O. We included Ga 3d electrons as valence states since it can be important for accurate determination of certain structural properties;^20,21 however, we found that it has minimal influence on the surface energies. To accurately describe the structural property and the electronic structure of monoclinic Al₂O₃ (denoted as θ-Al₂O₃) and β-Ga₂O₃, we use the hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE)²² with a mixing parameter 0.32, which produces lattice parameters and bandgaps in agreement with the experiment²⁰ (see the supplementary material). The structure of the bulk 10-atom monoclinic primitive cell of β-Ga₂O₃ and θ-Al₂O₃ is optimized with a 4 × 4 × 4 Monkhorst–Pack²³ k-point grid. The elastic constants of β-Ga₂O₃ and θ-Al₂O₃ are determined based on the conventional cell by performing six finite distortions of the lattice and deriving the elastic constants from the stress–strain relationship. Hellmann–Feynman forces for bulk calculations are converged to within 0.005 eV/Å. The elastic constants of the (Al_xGa_1−x)₂O₃ alloy are linearly interpolated between the end compounds. To quantitatively describe the bond strength, bond-stretching force constants are calculated using the finite difference method.

Surface calculations are performed in a slab geometry, with 12 Å vacuum separating the slabs. Details about the slabs used in the surface calculations are included in the supplementary material. For each slab, the top and bottom surfaces are equivalent by symmetry. For Brillouin-zone integrations, we used a 4 × 2 k-point grid in the plane and one k point in the out-of-plane direction. For unrelaxed surfaces, the atoms are frozen in their bulk positions. For relaxed surfaces, the Hellmann–Feynman forces are converged to within 0.01 eV/Å.

The conventional 20-atom unit cell of β-Ga₂O₃, with space group C2/m, is shown in Fig. 1(a). The correspondence between the conventional cell and the primitive cell was discussed in Ref. 24. Half of the Ga atoms are tetrahedrally coordinated Ga(I) and the other half are octahedrally coordinated Ga(II), as labeled in Fig. 1(a). Ga₂O₃ also has three types of O atoms: threefold coordinated O(I) (on a shared corner of two edge-sharing GaO₄ octahedra and one GaO₆ tetrahedron), threefold coordinated O(II) (on the shared corner of one GaO₆ octahedron and two GaO₄ tetrahedra), and fourfold coordinated O(III).

Figure 1 illustrates the six different surfaces considered here: (100)A, (100)B, (010), (001)A, (001)B, and (⁠$2¯$01). The surface orientations are defined with respect to the conventional cell. The labeling and definitions are consistent with those in Ref. 14. As will be discussed below, the unreconstructed surfaces are particularly stable because they satisfy electron counting. While surface reconstructions could occur on some surfaces under nonstoichiometric conditions during growth, they are not relevant for the crack-formation problem that is the focus of the current paper.

Atomic relaxations of the surface atoms are illustrated in Fig. 2. The magnitudes of atomic displacements drop dramatically for atoms away from the surface, showing that bulk-like positions are attained within a few atomic layers.

Table I lists the calculated surface energies for β-Ga₂O₃ and θ-Al₂O₃. Increasing the thickness of the slabs by including more atomic layers only changes the calculated surface energies by less than 0.01 J/m². Results from two previous calculations are included for comparison.^14,15 In general, our values agree reasonably well with those of Hinuma et al.;¹⁵ differences can be attributed to the fact that we use a more accurate functional. Differences with Ref. 14 are somewhat larger. Bermudez¹⁴ used a completely different computational approach (including functional, pseudopotentials, and basis set), which makes a detailed accounting for the differences difficult.

The surface energies are remarkably low—for instance, in comparison with GaN,²⁵ where even for the lowest-energy surface orientation [the (10$1¯$0) m plane], the surface energy is still 1.94 J/m². This high degree of stability for the β-Ga₂O₃ surfaces can be understood from the electron-counting rule.²⁶ For all surfaces considered here, the electrons transferred from dangling bonds on the Ga cations compensate the missing charges in the dangling bonds on the O anions, fulfilling the electron-counting rule and thereby lowering the surface energy.

We will be able to explain the relative stability of the various surfaces based on the following rules:

1. the nature of the broken bonds: softer bonds are easier to cut,

2. the dangling-bond density: lower densities of broken bonds reduce the surface energy, and

3. surface relaxation: atomic displacements can lower the energy by reducing strain and allowing rebonding.

The relative stability of the unrelaxed surfaces follows roughly the same trend in both β-Ga₂O₃ and θ-Al₂O₃. In β-Ga₂O₃, we have, in order of increasing surface energy, $Esurf(100)B<Esurf(100)A<Esurf(001)B⪅Esurf(010)⪅Esurf(2¯01)<Esurf(001)A$⁠. In θ-Al₂O₃, (001)B and (010) are interchanged, but their energies are close, anyway.

The most stable surface for both β-Ga₂O₃ and θ-Al₂O₃ is (100)B. This is consistent with the fact that the (100) plane is the easy cleavage plane in β-Ga₂O₃.²⁷ The low surface energy is mainly due to the fact that the (100)B surface cuts through easily broken Ga(II)–O(III) bonds in the GaO₆ octahedron (rule A). This can be quantified by inspecting the calculated bond-stretching force constants Φ_∥ listed in Table II. The table also provides information about the nature and density of bonds broken for each surface orientation. Table II shows that the bond-stretching force constant Φ_∥ for the Ga(II)–O(III) bond, with a value of 3.5 eV/Ų, is much smaller than that for any of the other bonds.

The (100)A surface also cuts Ga(II)–O bonds within a GaO₆ octahedron, producing the same dangling-bond density as the (100)B surface. However, while the (100)B surface cuts through soft Ga(II)–O(III) bonds, the (100)A surface cuts through stiffer Ga(II)–O(II) bonds (see Table II), resulting in a higher surface energy.

Similar to the (100)B surface, the (001)B surface also breaks soft Ga(II)–O(III) bonds. However, the dangling-bond density on the (001)B surface is twice as large as on (100)B (Table II), thus explaining the much higher surface energy (rule B).

On the other surfaces [i.e., (010), (⁠$2¯$01), and (001)A], bonds in a Ga(I)O₄ tetrahedron need to be broken. Since the Ga(I)–O bond in a GaO₄ tetrahedron is stronger than that in a Ga(II)O₆ octahedron (see Table II), this leads to higher surface energies (rule A).

When atomic relaxation is allowed, the relative stability of the relaxed surfaces follows the same trend in both β-Ga₂O₃ and θ-Al₂O₃. We have, in order of increasing surface energy, $Esurf(100)B<Esurf(100)A⪅Esurf(2¯01)<Esurf(001)B<Esurf(010)<Esurf(001)A$⁠. This trend agrees with the previous calculations,^14,15 although (as pointed out above) some numerical differences are evident.

Overall, the relative stability trend is very similar to that for unrelaxed surfaces, with one striking difference: the (⁠$2¯$01) surface is stabilized more than other surfaces, i.e., the relaxation energy for (⁠$2¯$01) is significantly larger than that for the other surfaces. This renders the energy of the relaxed (⁠$2¯$01) surface quite low, despite that bonds in Ga(I)O₄ tetrahedra need to be broken. Along the [$2¯$01] direction, the material is comprised of alternating layers of GaO₆ octahedra and GaO₄ tetrahedra [Fig. 1(c)]. The surface is formed by slicing through Ga–O bonds in the GaO₄ tetrahedra. Its unexpectedly low energy can be attributed to sizable atomic displacements, as visualized in Fig. 2(d); in fact, the relaxation energy is the largest among all of the investigated surfaces. The Ga(I) atom at the surface moves deeper into the layer, while the near-surface O(I) atom moves toward the surface. As a result, a new GaO₄ tetrahedron forms, with the original surface Ga(I) at its center; this tetrahedron is edge-sharing (instead of corner-sharing) with the neighboring GaO₆ octahedron. This dramatic structural relaxation (also noted in Ref. 15) greatly reduces the surface energy (rule C), as can be seen from the values in Table I.

Next, we calculate the elastic energy stored in the tensile-strained (Al_xGa_1−x)₂O₃ films based on E_elastic = 1/2∑_ijϵ_iC_ijϵ_j. For a monoclinic structure with the C2/m space group, there are 13 independent C_ij constants.²⁸ The calculated values are summarized in the supplementary material; they compare reasonably well with previous calculations²⁹ and with experiment.²⁹ 

We can now determine the elastic energy density in a thin film. First, we focus on pure θ-Al₂O₃ pseudomorphically grown on a β-Ga₂O₃ substrate (which could be either bulk Ga₂O₃ or a fully relaxed buffer layer). For a particular growth orientation, the in-plane components of the strain tensor are determined by the ratio of the in-plane lattice parameters of β-Ga₂O₃ and θ-Al₂O₃, while the out-of-plane components are obtained by minimizing the elastic energy. The resulting strains in θ-Al₂O₃ are listed in Table III for different growth orientations. The vanishing of the ϵ₄ and ϵ₆ components reflects the fact that for the growth scenarios discussed here, the unit cell remains monoclinic.¹⁰ 

For the growth of (Al_xGa_1−x)₂O₃ alloys, we follow the same procedure using linear interpolation between the end compounds to determine the in-plane components of ϵ and the elastic constants C_ij. Minimizing the elastic energy density at an Al fraction x then results in out-of-plane components of ϵ that are very nearly linear in x, indicating that they could equally be well directly linearly interpolated based on the strain tensor determined above for Al₂O₃ growth. As expected, the minimized elastic energy is quadratic in x.

We have presented all the information necessary to calculate the brittle fracture toughness and the elastic energy for (Al_xGa_1−x)₂O₃ films grown on a Ga₂O₃ substrate. We can now apply Eq. (1) to estimate a lower bound for h_c for each growth direction. The plane intersecting with the film surface that has the lowest surface energy is chosen to be the crack plane. The surface energy used in the Griffith theory should be the value for the unrelaxed surface, since the newly formed crack surfaces are infinitesimal portions of cleavage surfaces from the bulk crystal.³⁰ This affects the choice of crack plane. For example, for [100]-oriented film growth, both (⁠$2¯$01) and (001)B planes are possible crack planes. For relaxed surfaces, the energy of (⁠$2¯$01) is lower than that of (001)B (see Table I). However, we determine (001)B to be the proper crack plane, since its unrelaxed surface energy is lower than that of (⁠$2¯$01). For the [100], [010], and [001] film growth, the corresponding crack planes are (001)B, (100)B, and (100)B, respectively.

Figure 3 shows the calculated lower bound of critical thickness (h_c) as a function of Al concentration. Solid lines correspond to h_c using the surface energies for unrelaxed crack planes. Figure 3 shows that the [100] growth orientation gives rise to the largest h_c values for (Al_xGa_1−x)₂O₃ films; at x = 0.2, h_c is 157 nm. Critical thickness values for [010]- and [001]-oriented film growth are lower, yielding 34 nm and 27 nm at x = 0.2. The reduction of h_c for [010]- and [001]-oriented films is mainly due to the lower surface energy of the associated crack plane—(100)B, in both cases—in both β-Ga₂O₃ and θ-Al₂O₃.

We observe that critical layer thicknesses for (Al_xGa_1−x)₂O₃ films on β-Ga₂O₃ substrates are much smaller than those for Al_xGa_1−xN films on GaN substrates.³⁶ This is mainly due to the much lower surface energies for β-Ga₂O₃ and θ-Al₂O₃ compared with wurtzite GaN and AlN.³⁶ Additionally, while the elastic constants of β-Ga₂O₃ and θ-Al₂O₃ are comparable in magnitude to those of wurtzite GaN and AlN,³⁷ the lattice mismatch between the oxides is greater than that between the nitrides,³⁸ thus giving rise to larger elastic energies stored in (Al_xGa_1−x)₂O₃ films as resulting in smaller h_c.

To the best of our knowledge, for [010]-oriented uncracked (Al_xGa_1−x)₂O₃ film growth, the largest reported film thicknesses in samples grown by molecular beam epitaxy are from Ref. 31 and in samples grown by metal–organic chemical vapor deposition from Refs. 32–35 (see symbols in Fig. 3). While the calculated h_c at x = 0.4 is in excellent agreement with Ref. 35, the h_c at x ≤ 0.2 is underestimated compared to experiment.^31–33 We attribute this underestimation to the approximation made by replacing the fracture toughness Γ by brittle fracture toughness Γ_b in Eq. (1), i.e., the neglect of the energy of plastic deformation at the crack tip. This amounts to ignoring the atomistic nature of interactions at the crack tip. It has been shown based on a one-dimensional model³⁹ and atomistic simulations^30,40,41 that the discrete character of the lattice may lead to an increased force necessary to break a bond at the crack plane. This effect, known as “lattice trapping,” is neglected in our approximation using Γ_b. Lattice trapping contributes an additional energy barrier to propagate the crack and may lead to a higher critical thickness, particularly at lower Al concentrations where the elastic energy [denominator in Eq. (1)] is small.

In summary, we have calculated the surface energies of six surfaces of β-Ga₂O₃ and θ-Al₂O₃ and inferred the surface energies for (Al_xGa_1−x)₂O₃ alloys based on linear interpolation. We identified trends in the relative stability of various surfaces and explained these trends by considering (a) the nature of the broken bonds, (b) the dangling-bond density, and (c) the effect of relaxation. The (⁠$2¯$01) surface exhibits a pronounced atomic relaxation accompanied by a large reduction of the surface energy.

Employing the Griffith theory, the critical thickness for (Al_xGa_1−x)₂O₃ grown on Ga₂O₃ at different Al concentrations was calculated based on the elastic energy and the brittle fracture toughness. For (Al_xGa_1−x)₂O₃ films grown on Ga₂O₃ substrates, the crack planes identified to possess the lowest unrelaxed surface energy are found to be (001)B for [100] growth, (100)B for [010] growth, and (100)B for [001] growth. Our calculated values provide a lower bound to the actual critical thickness due to the neglect of lattice trapping (plastic deformation at the crack tip). Still, the numbers are in good agreement with the experiment. Our study provides insight into the physical mechanisms for crack formation in (Al_xGa_1−x)₂O₃ films grown on Ga₂O₃ substrates and the differences between various growth orientations.

See the supplementary material for details about structure of bulk β-Ga2O3 and θ-Al2O3, the slabs used in the surface calculations, and the calculated elastic constants.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors acknowledge fruitful discussions with Yuewei Zhang, James S. Speck, Hongping Zhao, Shuozhi Xu, and Mark Turiansky. This work was supported by the GAME MURI of the Air Force Office of Scientific Research (Grant No. FA9550-18-1-0479). Use was made of computational facilities purchased with funds from the National Science Foundation (NSF) (Grant No. CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; Grant No. NSF DMR 1720256) at UC Santa Barbara. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation under Grant No. ACI-1548562.