Crack formation limits the growth of (AlxGa1−x)2O3 epitaxial films on Ga2O3 substrates. We employ first-principles calculations to determine the brittle fracture toughness of such films for three growth orientations of the monoclinic structure: [100], [010], and [001]. Surface energies and elastic constants are computed for the end compounds—monoclinic Ga2O3 and Al2O3—and used to interpolate to (AlxGa1−x)2O3 alloys. The appropriate crack plane for each orientation is determined, and the corresponding critical thicknesses are calculated based on Griffith’s theory, which relies on the balance between elastic energy and surface energy. We obtain lower bounds for the critical thickness, which compare well with available experiments. We also perform an in-depth analysis of surface energies for both relaxed and unrelaxed surfaces, providing important insights into the factors that determine the relative stability of different surfaces. Our study provides physical insights into surface stability, crack planes, and the different degrees of crack formation in (AlxGa1−x)2O3 films for different growth orientations.

Due to its large bandgap (4.76 eV–5.04 eV),1–4 monoclinic β-Ga2O3 provides a large critical electric field strength in field effect transistors (FETs).5,6 Despite being a wide bandgap material, Ga2O3 can easily be n-type doped.5,7β-Ga2O3 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 (AlxGa1−x)2O3. For (AlxGa1−x)2O3 epilayers pseudomorphically grown on β-Ga2O3 substrates with the same crystallographic orientation as the substrate, tensile stress in (AlxGa1−x)2O3 is expected due to the lattice mismatch.10 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,11 V-shaped defects,11 or cracks,12 leading to a critical thickness where the strain can no longer be accommodated elastically. For the growth of β-(AlxGa1−x)2O3 films on Ga2O3 substrates, cracks have been found to be the major limitation.13 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 (AlxGa1−x)2O3 films grown on β-Ga2O3 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 surface14 is most stable and acts as the crack plane for (AlxGa1−x)2O3 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.

The Griffith criterion16 states that the critical thickness hc for a crack to propagate in an isotropic medium under equibiaxial stress is given by

hc=2ΓZ(i,jϵiCijϵj),
(1)

where Γ is the fracture toughness of the material (discussed below) and the sum in the denominator is over elastic constant (C) and strain (ϵ) components. Voigt notation is used, i.e., the components are defined as 1 = xx, 2 = yy, 3 = zz, 4 = yz, 5 = zx, and 6 = xy. Z is a dimensionless driving force that depends on the geometry of the crack and the film.12 In (AlxGa1−x)2O3/Ga2O3, “channeling” cracks have been observed,13 which initiate at the surface and extend to the substrate. For this channeling type, Z = 1.976.12 

Equation (1) represents an energy balance between the elastic energy stored in the strained film per unit area (h2i,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.12 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 hc constitutes a lower bound on the true critical thickness.

First-principles calculations were performed using DFT and projector augmented wave (PAW) potentials17 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 3d104s24p1 for Ga, 3s23p1 for Al, and 2s22p4 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 Al2O3 (denoted as θ-Al2O3) and β-Ga2O3, we use the hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE)22 with a mixing parameter 0.32, which produces lattice parameters and bandgaps in agreement with the experiment20 (see the supplementary material). The structure of the bulk 10-atom monoclinic primitive cell of β-Ga2O3 and θ-Al2O3 is optimized with a 4 × 4 × 4 Monkhorst–Pack23 k-point grid. The elastic constants of β-Ga2O3 and θ-Al2O3 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 (AlxGa1−x)2O3 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/Å.

Surface energies (Esurf) for stoichiometric slabs are calculated as

Esurf=12A(EslabNslabNbulkEbulk)=12Γb,
(2)

where A is the surface area, and Eslab and Ebulk are the total energies of the slab supercell and the bulk primitive cell. Nslab (Nbulk) is the number of formula unit (f.u.) contained in the slab (bulk) cells. The brittle fracture toughness (Γb) is twice the surface energy since two surfaces are formed, one on either side of the crack plane.16 Surface energies for (AlxGa1−x)2O3 are linearly interpolated between the end compounds.

The conventional 20-atom unit cell of β-Ga2O3, 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). Ga2O3 also has three types of O atoms: threefold coordinated O(I) (on a shared corner of two edge-sharing GaO4 octahedra and one GaO6 tetrahedron), threefold coordinated O(II) (on the shared corner of one GaO6 octahedron and two GaO4 tetrahedra), and fourfold coordinated O(III).

FIG. 1.

(a) The 20-atom conventional unit cell of monoclinic β-Ga2O3. The inequivalent Ga sites are labeled: tetrahedrally coordinated Ga(I) in green and octahedrally coordinated Ga(II) in purple and similarly for the O sites, threefold-coordinated O(I) in red, threefold coordinated O(II) in purple, and fourfold-coordinated O(III) in gray. The (100) and (001) surface cuts are depicted by black and blue dashed lines, respectively. “A” and “B” indicate different surface terminations. (b) The (010) surface cut (black dashed line). (c) The (2¯01) surface cut (black dashed line).

FIG. 1.

(a) The 20-atom conventional unit cell of monoclinic β-Ga2O3. The inequivalent Ga sites are labeled: tetrahedrally coordinated Ga(I) in green and octahedrally coordinated Ga(II) in purple and similarly for the O sites, threefold-coordinated O(I) in red, threefold coordinated O(II) in purple, and fourfold-coordinated O(III) in gray. The (100) and (001) surface cuts are depicted by black and blue dashed lines, respectively. “A” and “B” indicate different surface terminations. (b) The (010) surface cut (black dashed line). (c) The (2¯01) surface cut (black dashed line).

Close modal

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.

FIG. 2.

Atomic relaxations for (a) (100)B, (b) (010), (c) (001)B, and (d) (2¯01) surfaces of β-Ga2O3. The unrelaxed (left panel) and relaxed (right panel) surfaces are shown for each surface. Atomic displacements d relative to the ideal lattice sites are indicated in the left panels. For ease of visualization, the displacements are magnified by the factor listed in each panel. The blue circle in (d) depicts a newly formed tetrahedron after relaxation.

FIG. 2.

Atomic relaxations for (a) (100)B, (b) (010), (c) (001)B, and (d) (2¯01) surfaces of β-Ga2O3. The unrelaxed (left panel) and relaxed (right panel) surfaces are shown for each surface. Atomic displacements d relative to the ideal lattice sites are indicated in the left panels. For ease of visualization, the displacements are magnified by the factor listed in each panel. The blue circle in (d) depicts a newly formed tetrahedron after relaxation.

Close modal

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

TABLE I.

Surface energies Esurf (J/m2) for different surface orientations and terminations of β-Ga2O3 and θ-Al2O3 with and without atomic relaxations. Other calculations14,15 are listed for comparison.

β-Ga2O3θ-Al2O3
UnrelaxedRelaxedUnrelaxedRelaxed
SurfacesThis workReference 14 Reference 15 This workReference 14 Reference 15 This workReference 15 This workReference 15 
(100)B 0.60 0.96 0.61 0.34 0.68 0.38 1.06 0.74 0.62 0.45 
(100)A 1.39 1.68 1.28 0.85 1.13 0.80 2.07 1.68 1.26 1.06 
(010) 2.52 2.78 2.23 1.67 2.03 1.49 3.28 3.09 2.27 2.23 
(001)B 2.37 2.65 … 1.17 1.40 … 3.36 … 1.61 … 
(001)A 2.98 3.35 … 1.95 … … 3.78 … 2.56 … 
(2̄01) 2.67 … 2.16 0.96 … 0.75 3.52 2.94 1.28 1.04 
β-Ga2O3θ-Al2O3
UnrelaxedRelaxedUnrelaxedRelaxed
SurfacesThis workReference 14 Reference 15 This workReference 14 Reference 15 This workReference 15 This workReference 15 
(100)B 0.60 0.96 0.61 0.34 0.68 0.38 1.06 0.74 0.62 0.45 
(100)A 1.39 1.68 1.28 0.85 1.13 0.80 2.07 1.68 1.26 1.06 
(010) 2.52 2.78 2.23 1.67 2.03 1.49 3.28 3.09 2.27 2.23 
(001)B 2.37 2.65 … 1.17 1.40 … 3.36 … 1.61 … 
(001)A 2.98 3.35 … 1.95 … … 3.78 … 2.56 … 
(2̄01) 2.67 … 2.16 0.96 … 0.75 3.52 2.94 1.28 1.04 

The surface energies are remarkably low—for instance, in comparison with GaN,25 where even for the lowest-energy surface orientation [the (101¯0) m plane], the surface energy is still 1.94 J/m2. This high degree of stability for the β-Ga2O3 surfaces can be understood from the electron-counting rule.26 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.

1. Unrelaxed surfaces

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

The most stable surface for both β-Ga2O3 and θ-Al2O3 is (100)B. This is consistent with the fact that the (100) plane is the easy cleavage plane in β-Ga2O3.27 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 GaO6 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/Å2, is much smaller than that for any of the other bonds.

TABLE II.

Bond-stretching force constants (Φ, eV/Å2) in β-Ga2O3 between different cation–anion bonds. The number of bonds broken and the dangling-bond density (ρdb, Å−2) for different surfaces are also listed.

Φ(100)B(100)A(010)(2̄01)(001)B(001)A
Ga(I)–O(I) 13.4 
Ga(I)–O(II) 13.1 
Ga(I)–O(III) 10.5 
Ga(II)–O(I) 6.0 
Ga(II)–O(II) 5.9 
Ga(II)–O(III) 3.5 
ρdb  0.11 0.11 0.17 0.09 0.22 0.11 
Φ(100)B(100)A(010)(2̄01)(001)B(001)A
Ga(I)–O(I) 13.4 
Ga(I)–O(II) 13.1 
Ga(I)–O(III) 10.5 
Ga(II)–O(I) 6.0 
Ga(II)–O(II) 5.9 
Ga(II)–O(III) 3.5 
ρdb  0.11 0.11 0.17 0.09 0.22 0.11 

The (100)A surface also cuts Ga(II)–O bonds within a GaO6 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)O4 tetrahedron need to be broken. Since the Ga(I)–O bond in a GaO4 tetrahedron is stronger than that in a Ga(II)O6 octahedron (see Table II), this leads to higher surface energies (rule A).

2. Relaxed surfaces

When atomic relaxation is allowed, the relative stability of the relaxed surfaces follows the same trend in both β-Ga2O3 and θ-Al2O3. We have, in order of increasing surface energy, Esurf(100)B<Esurf(100)AEsurf(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)O4 tetrahedra need to be broken. Along the [2¯01] direction, the material is comprised of alternating layers of GaO6 octahedra and GaO4 tetrahedra [Fig. 1(c)]. The surface is formed by slicing through Ga–O bonds in the GaO4 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 GaO4 tetrahedron forms, with the original surface Ga(I) at its center; this tetrahedron is edge-sharing (instead of corner-sharing) with the neighboring GaO6 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 (AlxGa1−x)2O3 films based on Eelastic = 1/2∑ijϵiCijϵj. For a monoclinic structure with the C2/m space group, there are 13 independent Cij constants.28 The calculated values are summarized in the supplementary material; they compare reasonably well with previous calculations29 and with experiment.29 

We can now determine the elastic energy density in a thin film. First, we focus on pure θ-Al2O3 pseudomorphically grown on a β-Ga2O3 substrate (which could be either bulk Ga2O3 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 β-Ga2O3 and θ-Al2O3, while the out-of-plane components are obtained by minimizing the elastic energy. The resulting strains in θ-Al2O3 are listed in Table III for different growth orientations. The vanishing of the ϵ4 and ϵ6 components reflects the fact that for the growth scenarios discussed here, the unit cell remains monoclinic.10 

TABLE III.

Strain components for θ-Al2O3 grown on β-Ga2O3 substrates with [100], [010], and [001] orientations.

Film orientationϵ1 (%)ϵ2 (%)ϵ3 (%)ϵ4 (%)ϵ5 (%)ϵ6 (%)
[100] −3.8 4.5 3.6 −2.1 
[010] 3.7 −1.7 3.6 1.3 
[001] 3.7 4.5 −2.0 0.6 
Film orientationϵ1 (%)ϵ2 (%)ϵ3 (%)ϵ4 (%)ϵ5 (%)ϵ6 (%)
[100] −3.8 4.5 3.6 −2.1 
[010] 3.7 −1.7 3.6 1.3 
[001] 3.7 4.5 −2.0 0.6 

For the growth of (AlxGa1−x)2O3 alloys, we follow the same procedure using linear interpolation between the end compounds to determine the in-plane components of ϵ and the elastic constants Cij. 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 Al2O3 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 (AlxGa1−x)2O3 films grown on a Ga2O3 substrate. We can now apply Eq. (1) to estimate a lower bound for hc 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.30 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 (hc) as a function of Al concentration. Solid lines correspond to hc using the surface energies for unrelaxed crack planes. Figure 3 shows that the [100] growth orientation gives rise to the largest hc values for (AlxGa1−x)2O3 films; at x = 0.2, hc 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 hc 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 β-Ga2O3 and θ-Al2O3.

FIG. 3.

Calculated lower bound of the critical thickness (hc, nm) for (AlxGa1−x)2O3 grown on β-Ga2O3 as a function of Al concentration. (100) (blue), (010) (black), and (001) (red) orientated films are considered. Symbols denote the thickest uncracked films that have been reported for (010) growth: open black squares are from the study of Hilfiker et al.,31 the open black triangles and star are from the study of Bhuiyan et al.,32,33 the open circle is from the study of Ranga et al.,34 and the open diamond is from the study of Zhao.35 

FIG. 3.

Calculated lower bound of the critical thickness (hc, nm) for (AlxGa1−x)2O3 grown on β-Ga2O3 as a function of Al concentration. (100) (blue), (010) (black), and (001) (red) orientated films are considered. Symbols denote the thickest uncracked films that have been reported for (010) growth: open black squares are from the study of Hilfiker et al.,31 the open black triangles and star are from the study of Bhuiyan et al.,32,33 the open circle is from the study of Ranga et al.,34 and the open diamond is from the study of Zhao.35 

Close modal

We observe that critical layer thicknesses for (AlxGa1−x)2O3 films on β-Ga2O3 substrates are much smaller than those for AlxGa1−xN films on GaN substrates.36 This is mainly due to the much lower surface energies for β-Ga2O3 and θ-Al2O3 compared with wurtzite GaN and AlN.36 Additionally, while the elastic constants of β-Ga2O3 and θ-Al2O3 are comparable in magnitude to those of wurtzite GaN and AlN,37 the lattice mismatch between the oxides is greater than that between the nitrides,38 thus giving rise to larger elastic energies stored in (AlxGa1−x)2O3 films as resulting in smaller hc.

To the best of our knowledge, for [010]-oriented uncracked (AlxGa1−x)2O3 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 hc at x = 0.4 is in excellent agreement with Ref. 35, the hc 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 model39 and atomistic simulations30,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 β-Ga2O3 and θ-Al2O3 and inferred the surface energies for (AlxGa1−x)2O3 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 (AlxGa1−x)2O3 grown on Ga2O3 at different Al concentrations was calculated based on the elastic energy and the brittle fracture toughness. For (AlxGa1−x)2O3 films grown on Ga2O3 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 (AlxGa1−x)2O3 films grown on Ga2O3 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.

1.
H.
Tippins
,
Phys. Rev.
140
,
A316
(
1965
).
2.
T.
Matsumoto
,
M.
Aoki
,
A.
Kinoshita
, and
T.
Aono
,
Jpn. J. Appl. Phys., Part 1
13
,
1578
(
1974
).
3.
C.
Sturm
,
R.
Schmidt-Grund
,
C.
Kranert
,
J.
Furthmüller
,
F.
Bechstedt
, and
M.
Grundmann
,
Phys. Rev. B
94
,
035148
(
2016
).
4.
A.
Mock
,
R.
Korlacki
,
C.
Briley
,
V.
Darakchieva
,
B.
Monemar
,
Y.
Kumagai
,
K.
Goto
,
M.
Higashiwaki
, and
M.
Schubert
,
Phys. Rev. B
96
,
245205
(
2017
).
5.
M.
Higashiwaki
,
K.
Sasaki
,
A.
Kuramata
,
T.
Masui
, and
S.
Yamakoshi
,
Appl. Phys. Lett.
100
,
013504
(
2012
).
6.
W. S.
Hwang
,
A.
Verma
,
H.
Peelaers
,
V.
Protasenko
,
S.
Rouvimov
,
H.
Xing
,
A.
Seabaugh
,
W.
Haensch
,
C. G.
Van de Walle
,
Z.
Galazka
,
M.
Albrecht
,
R.
Fornari
, and
D.
Jena
,
Appl. Phys. Lett.
104
,
203111
(
2014
).
7.
J. B.
Varley
,
J. R.
Weber
,
A.
Janotti
, and
C. G.
Van de Walle
,
Appl. Phys. Lett.
97
,
142106
(
2010
).
8.
E.
Ahmadi
,
O. S.
Koksaldi
,
X.
Zheng
,
T.
Mates
,
Y.
Oshima
,
U. K.
Mishra
, and
J. S.
Speck
,
Appl. Phys. Express
10
,
071101
(
2017
).
9.
Y.
Zhang
,
A.
Neal
,
Z.
Xia
,
C.
Joishi
,
J. M.
Johnson
,
Y.
Zheng
,
S.
Bajaj
,
M.
Brenner
,
D.
Dorsey
,
K.
Chabak
,
G.
Jessen
,
J.
Hwang
,
S.
Mou
,
J. P.
Heremans
, and
S.
Rajan
,
Appl. Phys. Lett.
112
,
173502
(
2018
).
10.
M.
Grundmann
,
J. Appl. Phys.
124
,
185302
(
2018
).
11.
J.
Tersoff
and
F.
LeGoues
,
Phys. Rev. Lett.
72
,
3570
(
1994
).
12.
J.
Hutchison
and
T.
Wu
,
Advances in Applied Mechanics
(
Academic Press
,
San Diego, CA
,
1992
), Vol. 29.
13.
Y.
Zhang
and
J. S.
Speck
, private communication (
2020
).
14.
V.
Bermudez
,
Chem. Phys.
323
,
193
(
2006
).
15.
Y.
Hinuma
,
T.
Gake
, and
F.
Oba
,
Phys. Rev. Mater.
3
,
084605
(
2019
).
16.
A. A.
Griffith
,
Philos. Trans. R. Soc. London, Ser. A
221
,
163
(
1921
).
17.
P. E.
Blöchl
,
Phys. Rev. B
50
,
17953
(
1994
).
18.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
48
,
13115
(
1993
).
19.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
20.
H.
Peelaers
,
J. B.
Varley
,
J. S.
Speck
, and
C. G.
Van de Walle
,
Appl. Phys. Lett.
112
,
242101
(
2018
).
21.
S.
Mu
,
H.
Peelaers
, and
C. G.
Van de Walle
,
Appl. Phys. Lett.
115
,
242103
(
2019
).
22.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
118
,
8207
(
2003
).
23.
H. J.
Monkhorst
and
J. D.
Pack
,
Phys. Rev. B
13
,
5188
(
1976
).
24.
H.
Peelaers
and
C. G.
Van de Walle
,
Phys. Status Solidi B
252
,
828
(
2015
).
25.
C. E.
Dreyer
,
A.
Janotti
, and
C. G.
Van de Walle
,
Phys. Rev. B
89
,
081305
(
2014
).
26.
M.
Pashley
,
Phys. Rev. B
40
,
10481
(
1989
).
27.
Y.
Tomm
,
P.
Reiche
,
D.
Klimm
, and
T.
Fukuda
,
J. Cryst. Growth
220
,
510
(
2000
).
28.
S. V.
Gallego
,
J.
Etxebarria
,
L.
Elcoro
,
E. S.
Tasci
, and
J. M.
Perez-Mato
,
Acta Crystallogr., Sect. A: Found. Adv.
75
,
438
(
2019
).
29.
K.
Adachi
,
H.
Ogi
,
N.
Takeuchi
,
N.
Nakamura
,
H.
Watanabe
,
T.
Ito
, and
Y.
Ozaki
,
J. Appl. Phys.
124
,
085102
(
2018
).
30.
A.
Mattoni
,
L.
Colombo
, and
F.
Cleri
,
Phys. Rev. Lett.
95
,
115501
(
2005
).
31.
M.
Hilfiker
,
U.
Kilic
,
A.
Mock
,
V.
Darakchieva
,
S.
Knight
,
R.
Korlacki
,
A.
Mauze
,
Y.
Zhang
,
J.
Speck
, and
M.
Schubert
,
Appl. Phys. Lett.
114
,
231901
(
2019
).
32.
A. A. U.
Bhuiyan
,
Z.
Feng
,
J. M.
Johnson
,
Z.
Chen
,
H.-L.
Huang
,
J.
Hwang
, and
H.
Zhao
,
Appl. Phys. Lett.
115
,
120602
(
2019
).
33.
A. A. U.
Bhuiyan
,
Z.
Feng
,
J. M.
Johnson
,
H.-L.
Huang
,
J.
Sarker
,
M.
Zhu
,
M. R.
Karim
,
B.
Mazumder
,
J.
Hwang
, and
H.
Zhao
,
APL Mater.
8
,
031104
(
2020
).
34.
P.
Ranga
,
A.
Rishinaramangalam
,
J.
Varley
,
A.
Bhattacharyya
,
D.
Feezell
, and
S.
Krishnamoorthy
,
Appl. Phys. Express
12
,
111004
(
2019
).
35.
H. P.
Zhao
, private communication (
2020
).
36.
C. E.
Dreyer
,
A.
Janotti
, and
C. G.
Van de Walle
,
Appl. Phys. Lett.
106
,
212103
(
2015
).
37.
A.
Wright
,
J. Appl. Phys.
82
,
2833
(
1997
).
38.
H.
Schulz
and
K.
Thiemann
,
Solid State Commun.
23
,
815
(
1977
).
39.
R.
Thomson
,
C.
Hsieh
, and
V.
Rana
,
J. Appl. Phys.
42
,
3154
(
1971
).
40.
N.
Bernstein
and
D.
Hess
,
Phys. Rev. Lett.
91
,
025501
(
2003
).
41.
F.
Cleri
,
S. R.
Phillpot
,
D.
Wolf
, and
S.
Yip
,
J. Am. Ceram. Soc.
81
,
501
(
1998
).

Supplementary Material