Thermal conductivity of double polymorph Ga₂O₃ structures

Fabrication of heterostructures and functionalization of their properties is among the most successful strategies in solid-state technology, specifically for designing new electronic, magnetic, and thermal properties in semiconductors. There is a conventional way to fabricate such heterostructures by changing the chemical compositions of materials across interfaces, typically realized in situ during thin film synthesis.¹ Thus, there are numerous examples of highly useful heterostructures between chemically dissimilar semiconductors with tunable thermal properties.^2–5 On the other hand, the utilization of stacks between different phases in chemically identical materials, i.e., across homo-interfaces, is rare, perhaps with the exception of the crystalline/amorphous silicon structure exploited in high-efficiency solar cells.⁶ Nevertheless, the functional properties of the crystalline-to-crystalline polymorph junctions have not been investigated in a systematic way. This may be explained by the limited availability of such samples, which have both reasonable quality and technological relevance.

Indeed, there is practically no literature data reporting double polymorph structure synthesis with conventional deposition methods having technologically sufficient quality interfaces between polymorphs. This is quite understandable because the temperature/pressure growth conditions are strongly different for various polymorphs and must be changed abruptly, likely leading to interfacial quality degradation. Concurrently, even though polymorph structures may be induced by the localized pressure application to already existing crystals, control over the resulting structures is limited by the accuracy/scalability of the pressure application instrumentation.⁷ On the other hand, a new method to fabricate double γ-/β-Ga₂O₃ polymorph structures by a self-organized γ-polymorph transformation upon reaching a certain disorder threshold induced into β-Ga₂O₃ by irradiation was recently reported in the literature.^8–12 Importantly, the formation of such double polymorph γ-/β-Ga₂O₃ structures is characterized by the paradoxically good sharpness of the γ/β interface,¹⁰ making such structures comparable with conventional heterostructures in terms of functionalization options. As such, stacking different crystal symmetry polymorphs, in particular cubic γ-phase and monoclinic β-phase of Ga₂O₃, may result in step-like changes in thermal properties across these homo-interfaces, potentially exploitable in phononic devices in line with that of heterointerfaces between chemically dissimilar materials.^13–17 

Importantly, even though Ga₂O₃ attracts a lot of attention as a promising ultra-wide bandgap semiconductor, the understanding of its thermal transport properties is immature; e.g., experimental data are documented only for its stable β-Ga₂O₃ polymorph.^18–20 On the other hand, the thermal conductivity (k) data for the metastable Ga₂O₃ polymorphs are rare and limited to theoretical predictions, e.g., for the α-phase,^21–23, ε-phase,^22,24 and κ-phase.^23,25 This may be attributed to issues with the growth of stable polymorphs of reasonable quality. Thus, in this work, for the first time, we investigated the heat conduction properties of the double polymorph γ-/β-Ga₂O₃ structures, fabricated by a self-organized γ-polymorph transformation induced in β-Ga₂O₃ by irradiation.^9,10,12 For this reason, we used femtosecond laser-based time-domain thermoreflectance (TDTR) measurements with MHz modulation frequency to investigate heat conduction properties on a set of double γ/β polymorph Ga₂O₃ structures with variable top γ-layer thicknesses (ranging over hundreds of nanometers) fabricated on the β-Ga₂O₃ substrates having different crystallographic orientations. By comparing the experimentally measured data with the results of molecular dynamic (MD) simulations, we demonstrate up to an order of magnitude difference in thermal conductivities between γ- and β-phases depending on the crystalline orientation, providing additional options for thermal functionalization of the Ga₂O₃ devices.

As initial samples, we used high-purity, pristine β-Ga₂O₃ substrates with three different crystal orientations, such as (010), (−201), and (100), purchased from Novel Crystal Technology Inc. All samples had a lateral size of 5 × 5 mm² and were polished from the side exposed to the ion irradiation. In order to obtain γ-Ga₂O₃, all β-Ga₂O₃ samples were irradiated with 5 × 10¹⁶ Ga/cm² ions, which is known to be above the disorder threshold needed for converting the top part of the β-phase to the γ-phase at room temperature.¹² Notably, until recently, there was a puzzle with the identification of the new phase formed in the course of the disorder-induced ordering of β-Ga₂O₃. Some of the earlier works attributed the newly emerging polymorph to the orthorhombic κ-phase;^9,26 however, now it is unambiguously identified that the β-phase converts into the cubic spinel, i.e., the γ-phase.^8,10,12,27 Importantly, in order to vary the thickness of the top γ-phase layer, we used three different Ga ion energies, literally, 0.5, 1, and 1.7 MeV, at a constant fluence of 5 × 10¹⁶ Ga/cm².

To control the structural properties of the double Ga₂O₃ γ/β polymorph, we used a combination of x-ray diffraction (XRD) and Rutherford backscattering spectroscopy in channeling mode (RBS-C), relating this study to systematic structural characterization by extensive transmission electron microscopy imaging.^8–12,26,27 The XRD 2-theta scans were performed using a Bruker AXS D8 Discover diffractometer applying Cu Kα1 radiation in locked-coupled mode. The RBS-C analysis was carried out using 1.6 MeV He⁺ ions and was also used to measure the thickness of the newly formed top γ-layers. It is important to note that β-Ga₂O₃ is not amorphized under irradiation, but upon reaching a certain disorder threshold, it transforms into γ-phase, as shown previously,^9,10,12,26 with the characteristic box-like shape of the RBS-C spectra used as fingerprints of the double Ga₂O₃ γ/β polymorph in this study.

For the thermal transport measurements, we used a frequency-modulated (TDTR) setup; see details elsewhere.^28,29 Briefly, a Ti:sapphire mode-locked femtosecond laser (Tsunami, Spectra-Physics) at 782 nm wavelength, having an 80 MHz repetition rate and an 80 fs pulse duration, was used as the pump and the probe beams. The pump beam, modulated by an electro-optic modulator over 0.73–10 MHz, thermally excites the samples and controls heat penetration depth by $Dth=kπCvf$⁠, where k is the thermal conductivity, C_v is the volumetric specific heat capacity, and f is the pump modulation frequency. The probe beam is optically scan-delayed by a motorized delay stage and then detected by a photodetector connected to a radio-frequency lock-in amplifier. Both beams are focused on the sample surface by a 10× objective lens, resulting in a 1/e² diameter of ∼15 μm, which is substantially larger than the heat diffusion length. Therefore, TDTR measurements provide sensitivity only to cross-plane k at this beam diameter. A heat diffusion model was used to extract k and the interface thermal boundary conductance (G) of the studied samples by fitting a time-delayed signal ratio of the in-phase to out-of-phase voltages (–V_in/V_out) on a picosecond scale.²⁷ Notably, the Al heat transducing layer we deposited on all samples enabled the TDTR measurements. The thickness of this transducer was determined by the transient picosecond acoustic measurement and cross-checked with profilometer data, while the thermal conductivity of the Al transducer layer was measured using a reference sample; see the supplementary material.

To model k, we implemented classical molecular dynamics (MD) simulations using both equilibrium (EMD) and non-equilibrium (NEMD) molecular dynamics. For the EMD simulations, the lattice k was modeled with the Green–Kubo formalism relating the conductivity to the heat current autocorrelation function; see details in the supplementary material. NEMD modeling was also performed in comparison with our measured and EMD outcomes. We used reverse NEMD calculations,³⁰ where a heat source was set in the center of the simulation box with two heat sinks placed at opposing edges. The details are given in the supplementary material.

Figure 1 shows the summary of the experimental data for the phase-transformed Ga₂O₃, including (a) schematics of the sample preparation, (b)–(c) examples of the structural analysis, and (d) representative thermal conductivity data. Figure 1(b) demonstrates the XRD data of the samples fabricated on (010) oriented β-Ga₂O₃. There is a broad peak at ∼63.7°, in accordance with the literature identified as (440) cubic spinel γ-phase plane reflection forming as a result of the disorder-induced transition.^8–12,26,27 The increase in the γ-phase peak intensity along with the increase in the Ga⁺ ion energy is attributed to the broadening of the γ-phase film on top of the bulk β-Ga₂O₃ substrate, as illustrated in the schematics in Fig. 1(a) and as systematically confirmed in the literature.^8–12 For two other β-Ga₂O₃ crystal orientations, the evolutions of the β-to-γ-phase transition were similar to those in Figs. 1(a)–1(c) (not shown) and also consistently with earlier observations.^9,10,31

Figure 1(c) shows the RBS-C spectra of the (010)-oriented β-Ga₂O₃ samples irradiated with different energies, in comparison with the channeling and random data for the unimplanted sample. The top axis shows the corresponding depth scale. Importantly, we observe a tilted baseline for the unimplanted sample (as a signature of normal dechanneling) and the “box-like” RBS-C signal on the top of the tilted baseline for all irradiated samples as a result of the radiation disorder induced β-to-γ-phase transition.^9,10,31 The apparent broadening of the “box-like” part of the spectra as a function of the Ga⁺ ion energy is the evidence of the γ-phase layer thickening. The estimated γ-layer thickness in 0.5 MeV irradiated samples was ∼350 nm, while for 1 and 1.7 MeV, the γ-layer thickness proliferated to ∼650 and 1000 nm, respectively. This thickness variation is indicated by the vertical dashed lines in Fig. 1(c). Notably, none of the samples reached the “random” RBS-C level, confirming that β-to-γ-phase transformation is not accompanied by amorphization, consistently with the literature.^8–10,26,31 Heat conduction measurements also confirmed this RBS-C observation of not reaching the “amorphous limit,” as discussed below.

The greatest novelty of this work was measuring the thermal conductivity in the phase-transformed structures as illustrated in Fig. 1(d); however, before discussing these data, we have to introduce the TDTR measurements on pristine β-Ga₂O₃, for which we used a two-layer (Al/β-Ga₂O₃) model. The basic properties of β-Ga₂O₃, such as specific heat and density, were taken from the literature.^18–20 For reference, examples of the multi-frequency fitting to assess k values of β-Ga₂O₃ for different crystal orientations are available in the supplementary material. At this end, Fig. 1(d) plots the effective thermal conductivity (k_eff) as a function of f in the sample irradiated with 0.5 MeV, i.e., through the depth of the γ-/β-Ga₂O₃ structure. It is worth noting that k_eff was obtained assuming a two-layer model (Al/Ga₂O₃), which includes fitting the conductivity of the whole Ga₂O₃ layer and thermal interface conductance with a metal transducer layer. A non-uniform k_eff− trend is observed in Fig. 1(d) as a function of the heat penetration depth (D_th), nicely correlating with the double polymorph γ/β nature of the sample as illustrated by the sample cartoon in the inset in Fig. 1(d). Indeed, k_eff is the constant within the γ-layer and starts to increase at ∼1.9 MHz, corresponding to D_th ∼ 338 nm. Importantly, this depth is consistent with 350 nm γ-layer thickness as determined by RBS-C data in Fig. 1(c). Therefore, TDTR measurements also confirm the formation of a new layer with a lower thermal conductivity on the top of the β-Ga₂O₃ substrate. To assess k of γ-phase, we perform thermal data analysis based on a three-layer model (Al/γ-Ga₂O₃/β-Ga₂O₃). A similar analysis has been previously used for swift heavy ion irradiated single crystalline sapphire, where we used the same approach to spatially resolve the thermal conductivity of the radiation-induced subsurface amorphous Al₂O₃ layer and the ion track region of irradiated Al₂O₃.²⁸ Al and β-Ga₂O₃ properties were known from measurements of the unimplanted samples and were fixed during the fitting procedure for the samples containing γ-layers. Two unknown adjustable parameters, k of γ-Ga₂O₃ and conductance G between Al and γ-Ga₂O₃ phases, were obtained using the multi-frequency fitting (see the supplementary material). Importantly, in accordance with the literature, the interface between γ- and β-phases is remarkably abrupt.^8,9,27 Accounting that the sensitivity analysis showed minor variations of G between γ- and β-phases (see the supplementary material), in the current analysis G was set to a fixed value while the thickness of the γ/β interface was assumed to be ∼8 nm, as could be conservatively estimated from the literature.^9,12 The parameters used for fitting are summarized in Table SII in the supplementary material. Thus, for γ-phase, we determined k = 1.84 ± 0.1 W m⁻¹ K⁻¹ in the sample with the thinnest γ-film. This value is an order of magnitude lower than that in the bulk [010] oriented β-phase. However, it is larger than the thermal conductivity corresponding to the amorphous gallium oxide, estimated from Cahill’s minimum limit^32,33∼ 1.15 W m⁻¹ K⁻¹.

Figure 2 plots k_eff data as a function of the frequency for the variables (a) Ga⁺ ion energies and (b) crystal orientations. Importantly, in Fig. 2(a), we observe no k_eff variations in the samples irradiated by 1 and 1.7 MeV ions, as opposed to those in the 0.5 MeV irradiated sample. This is readily explainable by the γ-phase thickening to ∼650 and ∼1000 nm for samples irradiated with 1 and 1.7 MeV ions, respectively, as illustrated by the RBS-C data correlating with D_th not exceeding 500 nm at the lowest modulation rate of 0.73 MHz. It implies that with TDTR at this modulation frequency, we measure k_eff exclusively of the γ-phase but not of the underlying β-phase. The results of the multi-layer fitting applied to the 350, 650, and 1000 nm thick γ-films, corresponding to the irradiations with 0.5, 1, and 1.7 MeV, ions respectively, are summarized in Table I. Notably, for all previously studied semiconductors, one would expect a drop in k with the increase in ion energy because of increased defect production. However, the results in Fig. 2(a) show a rather constant trend, perhaps exhibiting a slight increase in k as a function of the irradiation energy, which can be attributed to even better stabilization of the γ-phase obtained via disorder-induced ordering in the higher energy irradiated samples.

Furthermore, Fig. 2(b) confirms that k of the newly formed γ-phase is independent of the crystal orientation and remains at ∼ 1.8 W m⁻¹ K⁻¹. These results indicate that the phase transformation in β-Ga₂O₃ occurs independent of the crystal orientation, consistently with the literature.^9,31 Another important correlation is that k of the cubic spinel γ-phase is isotropic, in contrast to the β-phase anisotropy, as also illustrated in Fig. 2(b) by comparing k values for different crystallographic directions in β-Ga₂O₃.

Furthermore, using equilibrium molecular dynamics (EMD) simulations, we calculated the k values of β-Ga₂O₃ for three different crystal directions; see the supplementary material. Notably, significant oscillations in the autocorrelation function while using the Green–Kubo relation inherently contribute to high uncertainties for k because of the high noise levels preventing the identification of the convergence region. Therefore, 20 independent trajectories were chosen to improve the statistics so that the conductivity reached equilibrium after 240 ps of correlation time. These calculated β-Ga₂O₃ thermal conductivity results using the Born–Mayer–Huggins (BMH) potential³⁴ are in good agreement with the literature data^18,19,34 exhibiting a strong anisotropy with the largest conductivity along the [010] direction (see the supplementary material). On the other hand, the tabGAP potential³⁵ underestimates the k value for β-phase, resulting in 11.3 ± 0.7 Wm⁻¹K⁻¹ along the [010] direction, which is almost two times lower than the Born–Mayer–Huggins potential and experimental values reported before (see the supplementary material). Nevertheless, using the same approach and two different potentials, we performed calculations of the k- in the γ-phase; see the data in Fig. 3(a). For these simulations, the saturation of the k-values was achieved at 160 ps, resulting in k of 4 W m⁻¹ K⁻¹ for the Born–Mayer–Huggins potential and 3.1 W m⁻¹ K⁻¹ for the tabGAP potential, i.e., showing a qualitatively consistent trend with the experimental values in Fig. 2 and Table I; however, exceeding the experimental values approximately by a factor of two.

For more insights into mechanisms, we performed complementary NEMD simulations using Born–Mayer–Huggins potential; in particular, Fig. 3(b) shows the evolution of 1/k as a function of the inverse size of the simulation box in γ-Ga₂O₃. Figure 3(b) shows an apparent linear trend, with the intercept of the 1/k axis giving the k value for an “infinitely large” box corresponding to a bulk thermal conductivity of 3.93 W m⁻¹ K⁻¹ for γ-Ga₂O₃ consistently with EMD results.

At this end, Fig. 4(a) summarizes the k-values for β- and γ-Ga₂O₃ obtained in the present work in comparison with the literature data for different forms of Ga₂O₃.^{21–25,33,36,37} Most importantly, our experimental and simulation results show that the thermal conductivity of γ-phase is significantly lower than that of other Ga₂O₃ crystal phases, which can be attributed to a lower phonon mean free path λ_mfp ≈ 4 nm in γ-phase as extracted from the NEMD simulations. Notably, the NEMD results show that the k-values of the γ-phase are quite independent of the film thickness in comparison with those in β-Ga₂O₃ thin films.³⁶ This is comparable with the trend observed in β-(Al_xGa_1−x)₂O₃ thin films with different Al contents, where the k variations were mainly attributed to the phonon-alloy disorder scattering,³⁷ implying that the thermal transport in these structures is dominated mostly by the vibrational modes with short λ_mfp. Notably, it was reported that β-(Al_xGa_1−x)₂O₃ converts to its γ-phase when the Al concentration exceeds 40%.^38,39

Importantly, as already mentioned above, our experimentally and theoretically assessed k-values in γ-Ga₂O₃ are in qualitative agreement; however, the experimental data in Fig. 2 are approximately twofold lower than the MD data. This difference may be attributed to the fact that the MD simulations in Fig. 3 have not accounted for a potential strain accumulation in the double γ-/β-Ga₂O₃ samples, as was discussed in the literature.⁹ Thus, to check whether the strain effect may explain the discrepancy, we repeated the EMD simulations, applying in-plane (ε_xx) and out-of-plane (ε_yy) strain in ratios similar to those reported in the literature.⁹ The corresponding simulation calculations show, see the supplementary material, that the k-values of the γ-phase exhibit a decreasing trend as a function of the applied strain reaching the experimentally measured k value at ε_xx = 6.3% and ε_yy = 1.86%. Notably, even though the impact of strain, as discussed above, seems reasonable in the context of the literature, at the present stage, we cannot rule out other reasons for the reduction of heat conductivity, e.g., associated with ion-induced point defects and/or the imperfections in the γ-films.⁴⁰ Nevertheless, it also should be noted that the double γ/β polymorph structures exhibit remarkably high radiation tolerance for consequent radiations,^10,41 i.e., the residual defect level in γ-Ga₂O₃ film is comparatively low compared with conventional cases of radiation defect accumulation leading to amorphization in materials.

This work paves a new direction for studying nanoscale heat propagation in double γ-/β-Ga₂O₃ polymorph structures, potentially applicable to modulate functional thermal properties across interfaces without changing the chemical composition of materials, in contrast to that in conventional heterostructures. Consistently with the literature, the samples were fabricated by a self-organized β-to-γ polymorph transformation upon reaching a certain disorder threshold induced into β-Ga₂O₃ by the ion irradiation, as controlled by a combination of the XRD and RBS-C. To make a systematic comparison of the thermal properties, we depth-profiled the thermal conductivity across γ-/β-Ga₂O₃ structures having variable γ-layer thicknesses (350–1000 nm) on top of the β-Ga₂O₃ substrates using TDTR. As a result, in γ-Ga₂O₃, k was in the range between 1.84 and 2.11 W m⁻¹ K⁻¹ for 350 and 1000 nm samples, respectively, and independently of the initial β-substrates orientations. In parallel, we performed heat conduction MD simulations in the γ- and β-Ga₂O₃ lattices, using both equilibrium and non-equilibrium approaches. Importantly, for the relaxed γ-Ga₂O₃ lattice, in the MD simulations, we obtained approximately twofold higher k in comparison with the experimental results. Consistently with the literature, such discrepancy may be attributed to the residual biaxial strain remaining in the γ-Ga₂O₃ lattice after the β-to-γ transition, even though other reasons for the reduction of thermal conductivity associated with the imperfections in the new γ-films could not be ruled out. In its turn, for monoclinic β-Ga₂O₃, the heat propagation showed a distinct anisotropy, with values ranging from 10 W m⁻¹ K⁻¹ for [−201] to 20 W m⁻¹ K⁻¹ for [010] orientations, consistently with the MD simulations. Thus, these results demonstrate the viability of the nanoscale thermoreflectance depth-profiling metrology across chemically identical polymorph interfaces and open an avenue to explore the nanoscale phonon transport in such “polymorph heterostructures” exhibiting variations in thermal conductivities, e.g., of an order of magnitude in γ-/β-Ga₂O₃ formed on [010] oriented β-substrates.

The details of data analysis and molecular dynamics simulations are given in the supplementary material.

This work was supported by the Science Committee of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant Nos. AP19577063 and AP19679332), Nazarbayev University grants via the Collaborative Research Program (CRP) (Grant No. 11022021CRP1504), the Faculty Development Competitive Research Grants Program (FDCRGP) (Grant No. 20122022FD4130), and the M-ERA.NET GOFIB project (Research Council of Norway Project No. 337627). The international collaboration was in part enabled by the INTPART Program funded by the Research Council of Norway (Project No. 322382).

The authors have no conflicts to disclose.

Azat Abdullaev: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Kairolla Sekerbayev: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal). Alexander Azarov: Formal analysis (equal); Investigation (equal); Validation (equal). Vishnukanthan Venkatachalapathy: Formal analysis (supporting); Investigation (supporting); Validation (supporting). Vinay S. Chauhan: Formal analysis (equal); Investigation (equal); Validation (equal). Zhandos Utegulov: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Resources (equal); Writing – review & editing (equal). Andrej Kuznetsov: Conceptualization (equal); Methodology (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.