# Spectral Thermal Spreading Resistance of Wide-Bandgap Semiconductors in Ballistic-Diffusive Regime

Yang Shen<sup>D</sup>, Yu-Chao Hua, Han-Ling Li<sup>D</sup>, S. L. Sobolev<sup>D</sup>, and Bing-Yang Cao<sup>D</sup>

**Abstract**— To develop efficient thermal management strategies for wide-bandgap (WBG) semiconductor devices, it is essential to have a clear understanding of the heat transport process within the device and accurately predict the junction temperature. In this article, we use the phonon Monte Carlo (MC) method with the phonon dispersion of several typical WBG semiconductors, including GaN, SiC, AlN, and β-Ga₂O₃, to investigate the thermal spreading resistance in a ballistic-diffusive regime. This work shows that when compared with Fourier's law-based predictions, the increase in the thermal resistance caused by the ballistic effects is strongly related to the phonon dispersion. Based on the model derived under the gray-medium approximation and the results of dispersion MC, we obtained a thermal resistance model that can well address the issues of thermal spreading, ballistic effects, and the influence of phonon dispersion. The model can be easily coupled with finite-element method (FEM)-based thermal analysis and applied to different materials. This article can provide a clearer understanding of the influence of phonon dispersion on the thermal transport process, and it can be useful for the prediction of junction temperatures and the development of thermal management strategies for WBG semiconductor devices.

**Index Terms**— Ballistic transport, phonon Monte Carlo (MC) simulation, thermal spreading resistance, wide-bandgap (WBG) semiconductor.

Manuscript received January 23, 2022; revised March 16, 2022; accepted April 12, 2022. Date of publication April 29, 2022; date of current version May 24, 2022. This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grant 51825601, Grant U20A20301, and Grant 20-58-53017; and in part by the Russian Foundation for Basic Research (RFBR) under Grant 52011530030. The review of this article was arranged by Editor S.-M. Hong. (Corresponding author: Bing-Yang Cao.)

Yang Shen, Han-Ling Li, and Bing-Yang Cao are with the Key Laboratory of Thermal Science and Power Engineering of Education of Ministry, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China (e-mail: sy980829@163.com; lihanling1994@163.com; caoby@tsinghua.edu.cn).

Yu-Chao Hua is with the LTEN Laboratory, UMR 6607, Polytech Nantes, University of Nantes, F-44000 Nantes, France (e-mail: huayuchao19@163.com).

S. L. Sobolev is with the Institute of Problems of Chemical Physics, Academy of Sciences of Russia, 142432 Chernogolovka, Russia, and also with the Heat Power Department, Samara State Technical University, 443100 Samara, Russia (e-mail: sobolev@icp.ac.ru).

Color versions of one or more figures in this article are available at
https://doi.org/10.1109/TED.2022.3168798

Digital Object Identifier 10.1109/TED.2022.3168798

## I. INTRODUCTION

WIDE-BANDGAP (WBG) semiconductor devices such as GaN high-electron-mobility transistors (HEMTs) [1], [2], SiC metal-semiconductor field-effect transistors (MESFETs) [3], [4], and recently burgeoning AlN and β-Ga₂O₃-based ultra-WBG (UWBG) devices [5]–[7] are extremely attractive for high-power and high-frequency electronic applications, due to the high breakdown voltage led by the WBG. However, owing to their high power density, these devices usually hold very high junction temperatures. The significant overheating within the devices largely restricts their actual performance [8]–[10] and shortens the device lifetime [11]. Therefore, it is important to have a clear understanding of the heat transport process within the device, to accurately predict the junction temperature and therefore develop effective thermal management strategies [12]–[14].

In general, there are two prominent features of the heat transport process within WBG semiconductor devices. First, the heat generation is concentrated at the top of the channel layer [15], [16]. The heat source region is extremely small [17] compared with the channel layer's length and width. When heat spreads from a small heat source region to a much larger area, there is a significant thermal spreading resistance, which is one of the primary sources of the total thermal resistance [18] and finally results in a large near-junction temperature spike [19]. Second, the basic structure of the device is made up of multilayer micro/nanofilms, and the thickness of the channel layer is usually less than 3 µm [16], [20], [21]. This characteristic size is comparable with the mean free paths (MFPs) of phonons, which are the main heat carriers in semiconductors [22]. The transport of phonons in the device is significantly suppressed by strong phonon-boundary scatterings, leading to an additional ballistic thermal resistance [23] and further increasing the junction temperature [24]. In this case, Fourier's law becomes inapplicable and ballistic-diffusive heat conduction emerges [25]. Therefore, to accurately predict the junction temperature, it is important to provide a clear understanding of the thermal spreading process and the influence of ballistic effects.

Although great efforts have been devoted to multiscale simulation techniques to predict the temperature profile of the device [26]–[28], little work has been done for detailed analysis of the thermal spreading process and the influence
of ballistic effects [24]. In addition, these simulations often take a long time to converge, and when considering phonon dispersion, they can be much more time-consuming [29]. The high computational demands partly restrict the applications of these simulation techniques in the practical 3-D thermal analysis and optimizations of real devices, where the finite-element method (FEM) still plays a major role [30]–[32]. Therefore, it is still urgent and significant to propose new physical insights into the thermal transport process within the device and therefore develop computationally efficient models that can take multiple factors into account and can give a fast estimation of the junction temperature. Additionally, the model should be easily coupled with the widely used FEM-based design process.

Based on Fourier's heat conduction law, analytical thermal spreading resistance models have been extensively studied [18], [19], [33]. Recently, Hua et al. [24] developed a semiempirical thermal resistance model that could consider the influences of both thermal spreading and ballistic effects under the gray-medium approximation. The work indicates that ballistic effects can significantly alter the temperature distributions within the channel layer and can greatly increase the thermal resistance. However, this work did not scrutinize the influence of phonon dispersion. Since thermal transport in semiconductors is the contribution of all phonon modes, different phonon dispersion relations will significantly affect the junction temperature [29]. For typical WBG semiconductors, Freedman et al. [34] found that phonons with MFPs greater than 1000 ± 230 nm, 2500 ± 800 nm, and 4200 ± 850 nm contribute 50% of the bulk thermal conductivity of GaN, AlN, and 4H-SiC near room temperature, respectively. In consideration of the characterized size of the channel layer, the phonon transport in the device should undergo strong size effects caused by the wide distributions of phonon MFPs. Therefore, it is necessary to reexamine the applicability of the semiempirical gray model and analyze the impacts of phonon dispersion in more detail.

In this work, we use phonon Monte Carlo (MC) methods with the phonon dispersion of several typical WBG semiconductors to study the thermal spreading resistance in a ballistic-diffusive regime. The considered materials include GaN, SiC, AlN, and β-Ga₂O₃. It was found that the gray model was insufficient to reflect the strong ballistic effects caused by the wide phonon MFP distributions. Based on the analysis of the deviations, material-independent correction factors were introduced to the gray model to better reflect the influence of phonon dispersion on the thermal transport process. The improved model can well address the issues of thermal spreading, ballistic effects, and the influence of phonon dispersion, and it is easy to be coupled with FEM-based thermal analysis. Our work could be helpful for junction temperature predictions and the development of thermal management strategies for WBG semiconductors.

## II. PROBLEM FORMULATION AND METHODOLOGIES

### A. Problem Statement

Fig. 1 shows the basic system for thermal spreading resistance analysis, which is a representative 2-D simplified model

<figure><img src="image_1.png" /><figcaption>Fig. 1. Schematic for the basic system for thermal spreading resistance analysis.</figcaption></figure>

for the channel layer. The geometry of the system can be represented by three parameters: the channel thickness $t$, the width of the structure $w$, and the width of the heating area $w_g$. At the top surface of the structure, uniform heat flux is specified over the region of the heat source, and the total heating power is denoted by $Q$. The remainder of the heat source surface is taken as adiabatic. The lateral boundaries are set as periodic, and the bottom boundary serves as an isothermal heat sink.

The total thermal resistance can be calculated by using the mean temperature over the heat source area ($\bar{T}_s$), heat sink temperature ($T_0$), and total heating power ($Q$) [18]

$$R_t = \frac{\bar{T}_s - T_0}{Q} \qquad (1)$$

For convenience in identifying the thermal resistance model [24], the dimensionless total thermal resistance is defined as $R_t/R_{1d_0}$. $R_{1d_0}$ is the 1-D thermal resistance of purely diffusive heat conduction

$$R_{1d_0} = t/(wk_0) \qquad (2)$$

where $k_0$ is the intrinsic thermal conductivity.

### B. Semiempirical Thermal Resistance Model

For the system shown in Fig. 1, the total thermal resistance can be attributed to three parts: the thermal spreading process, the cross-plane phonon ballistic transport, and the ballistic effect with $w_g$ comparable with MFP [24]. Thus, the total thermal resistance can be rearranged to

$$\frac{R_t}{R_{1d_0}} = \frac{R_F}{R_{1d_0}} \cdot \frac{R_{1d}}{R_{1d_0}} \cdot \left[ \frac{R_t}{R_{1d}} \left( \frac{R_F}{R_{1d_0}} \right)^{-1} \right] \qquad (3)$$

where

$$\frac{R_F}{R_{1d_0}} = 1 + \left(\frac{w}{w_g}\right)^2 \left(\frac{w}{t}\right)^\infty \sum_{n=1}^\infty \frac{8 \sin^2\left(\frac{w_g n \pi}{2w}\right) \cdot \cos^2\left(\frac{n \pi}{2}\right)}{(n \pi)^3 \coth\left(\frac{n \pi}{w}\right)} \qquad (4)$$

$$\frac{R_{1d}}{R_{1d_0}} = 1 + \frac{2}{3}Kn_t \qquad (5)$$

$$\frac{R_t}{R_{1d}} \left( \frac{R_F}{R_{1d_0}} \right)^{-1} = r_w = 1 + A_w (w_g/w, w/t) K n_w \qquad (6)$$

$R_F/R_{1d_0}$ is the analytical expression of the dimensionless total thermal resistance derived based on Fourier's law [35], which is used to characterize the thermal spreading effect.
The remaining two terms are used to reflect the ballistic thermal resistance. To characterize the strength of the ballistic effects, two Knudsen numbers $Kn_t$ and $Kn_w$ are defined as

$$Kn_t = l_0/t \qquad (7)$$

$$Kn_w = l_0/w_g \qquad (8)$$

where $l_0$ is the intrinsic phonon MFP. $R_{1d}/R_{1d,0}$ can be derived with the differential approximation and the temperature jump boundary conditions [36], which represents the cross-plane ballistic thermal resistance. The last term denoted by $r_w$ is used to reflect the ballistic effect with $w_g$ comparable with MFP. $A_w$ is a fitting parameter that is a function of geometric factors $w_g/w$ and $w/t$.

The model is derived under the gray-medium approximation, that is, for the single-phonon mode. When considering phonon dispersion, the size effects on different phonon modes are different. In this case, the effective thermal conductivity can be evaluated by using an integral with the mode-dependent modified MFPs

$$k_{\text{eff}} = \frac{1}{3} \sum_{j} \int_{0}^{\omega_{j}} \hbar \omega \frac{\partial f_{0}}{\partial T} \text{DOS}_{j}(\omega) v_{g, \omega, j} l_{m, j} d\omega \quad (9)$$

where

$$l_{m,j} = \frac{l_{0,j}}{(1 + A_w Kn_{w-\omega,j})(1 + \frac{2}{3} Kn_{t-\omega,j})} \quad (10)$$

in which $l_{0,j}$ is the intrinsic frequency-dependent MFP of the j phonon branch, and $Kn_{w-\omega,j} = l_{0,j}/w_g$ and $Kn_{t-\omega,j} = l_{0,j}/t$ are the frequency-dependent and phonon branch-dependent Knudsen numbers, respectively. With the effective thermal conductivity, the total thermal resistance can be calculated with Fourier's law-based model as (4).

### C. Phonon Monte Carlo Simulations

In this article, the phonon-tracing MC technique is used to simulate ballistic-diffusive heat conduction [37], [38]. The basic simulation settings are the same as those in [24]. The main difference is that to consider the phonon dispersion, the emitted phonon bundles are sampled from the phonon spectrum, and their properties are redetermined after phonon-phonon scattering [29]. In this work, the energy-based variance-reduced technique proposed by [39] is adopted to determine the probability of drawing a phonon bundle at a certain polarization and frequency. Details of dispersion MC can be found in [29].

Four typical WBG or UWBG semiconductors are considered in this work, including GaN, AlN, SiC, and β-Ga₂O₃. An isotropic sine-shaped phonon dispersion (Born-von Karman dispersion) is used for all materials [40], $\omega(k) = \omega_{\text{max}} \sin(\pi k/2 k_m)$, where $k_m = (6\pi^2 n)^{1/3}$ with n as the volumetric density of primitive cells. The previous research has verified that this model can well reflect the phonon MFP spectrum of different materials [34]. Relaxation time is also needed by dispersion MC to determine the MFPs of different phonon modes. The essential phonon scattering mechanisms of these materials include impurity scattering (I) and Umklapp phonon-phonon scattering (U). The relaxation time can be expressed as $\tau_I^{-1} = A\omega^4$ and $\tau_U^{-1} = B\omega^2T \exp(-C/T)$ [41],

<table><caption>TABLE I
PHONON DISPERSION AND RELAXATION TIME PARAMETERS FOR DIFFERENT SEMICONDUCTOR MATERIALS</caption><thead><tr><th>Parameter (Unit)</th><th>GaN</th><th>AlN</th><th>SiC</th><th>β-Ga₂O₃</th></tr></thead><tbody><tr><td>$k_0$ (1 × 10⁹ m⁻¹)</td><td>10.94</td><td>11.19</td><td>8.94</td><td>6.74</td></tr><tr><td>$\omega_m$ (1 × 10¹³ rad/s)</td><td>3.50</td><td>5.18</td><td>7.12</td><td>1.6</td></tr><tr><td>$a_D$ (Å)</td><td>2.87</td><td>2.81</td><td>3.51</td><td>4.66</td></tr><tr><td>$A$ (1 × 10⁻⁴⁵ s³)</td><td>5.26</td><td>10.5</td><td>1.00</td><td>1.38E-6</td></tr><tr><td>$B$ (1 × 10⁻¹⁹ s/K)</td><td>1.10</td><td>0.728</td><td>0.596</td><td>9.31</td></tr><tr><td>$C$ (K)</td><td>200</td><td>287.5</td><td>235.0</td><td>62.6</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Fig. 2. Thermal conductivity from the model calculations (line), and from experiments (symbols) of GaN [45], AlN [46], 6H-SiC [47], and β-Ga₂O₃ [42], [44].</figcaption></figure>

where A, B, and C are the fitting constants. The total relaxation time can be calculated using Matthiessen's rule $\tau^{-1} = \tau_I^{-1} + \tau_U^{-1}$.

Parameters of different materials can be obtained by fitting the measured thermal conductivities. Table I lists the parameters used in dispersion MC in this work. For GaN, AlN, and SiC, the fitted values in [27] are used. For β-Ga₂O₃, the parameters are obtained by fitting the measured thermal conductivities along the [100] crystallographic direction [42]–[44]. As shown in Fig. 2, the model can well reflect the variance of the thermal conductivity with the temperature of different materials.

## III. RESULTS AND DISCUSSION

### A. Total Thermal Resistance

Fig. 3 compares the values of dimensionless total thermal resistance of GaN calculated by dispersion MC and the semiempirical thermal resistance model. The predictions of Fourier's law-based model are also given for comparison. The results with different $w_g/w$ and of other materials are not exhibited since they show similar patterns. From Fig. 3, it can be seen that both the model predictions and MC simulation results are much higher than the results predicted based on Fourier's law, which indicates that the ballistic effects significantly increase the thermal resistance.

However, it can be noted that there are nonnegligible deviations between the dimensionless total thermal resistance predicted by the model and dispersion MC. For GaN, in particular, the maximum relative deviation achieves 27% with $t = 4$ µm and exceeds 33% with $t = 0.2$ µm. As shown
<figure><img src="image_3.png" /><figcaption>Fig. 3. Dimensionless total thermal resistance of GaN as a function of w/t predicted by the model and MC simulations, with $w_g/w = 0.01$, $t = 0.2$, and 4 µm, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4. Dimensionless total thermal resistance of GaN as a function of w/t predicted by the model based on $l_{m,j,1}$ and MC simulations, with $w_g/w = 0.01$, $t = 0.2$ and 4 µm, respectively.</figcaption></figure>

in Fig. 3, the deviations can be divided into two categories. When $w/t$ is small, the model predicted thermal resistance tends to be larger than those of dispersion MC. Moreover, with increasing $w/t$, both the dimensionless thermal resistance predicted by the model and dispersion MC increase and reach a plateau value. In the plateau region, the model predictions are lower than the MC simulation results.

The deviations indicate that the gray model is insufficient to reflect the wide span of the phonon spectrum. To figure out the source of the deviations, the key issue is to separate the impacts of the different ballistic effects. In the following, we first focus on the deviations related to the cross-plane ballistic part. Then, by offsetting the discrepancy caused by the cross-plane ballistic effect, the $Kn_w$ dependence of the deviations from the ballistic effect with $w_g$ comparable with MFP is illustrated.

## B. Cross-Plane Ballistic Effect

From (10), it can be found that in the gray model, different ballistic effects are embodied in different terms of the modified MFPs. By substituting $(1+A_w(w_g/w, w/t)Kn_{w-\omega,j})$ with the $r_w$ directly calculated from dispersion MC, as shown in the following equation:

$$l_{m,j,1} = \frac{l_{0,j}}{(1 + \frac{2}{3}Kn_{t-\omega,j})r_{w,\text{dispersion}}} \quad (11)$$

and calculating the effective thermal conductivity and the total thermal resistance using $l_{m,j,1}$, the differences between model predictions and MC simulation results caused by the ballistic effect with $w_g$ comparable with MFP can be compensated. Thus, the remaining deviations should be attributed to the cross-plane ballistic effect.

Fig. 4 shows the dimensionless total thermal resistance predicted by dispersion MC and the model based on $l_{j,m,1}$. In this case, the model predictions are always lower than the MC simulation results. Comparing Fig. 4 with Fig. 3, it can be found that the anomalous tendency of the model predicted thermal resistance with small $w/t$ disappears. This indicates that the deviations between the thermal resistance

predicted by the original model and dispersion MC with small $w/t$ are mainly from the ballistic effect with $w_g$ comparable with MFP. In addition, as shown in Fig. 4, the deviations increase significantly as $t$ decreases from 4 µm to 0.2 µm. This suggests that the gray model's capability to reflect the cross-plane ballistic effect is further impaired with decreasing channel thickness. To depict the channel thickness dependence of the model's insufficiency, an appropriate average MFP must be defined to represent the entire phonon spectrum of different materials.

Li et al. [29] compared different approaches to calculate the average MFP. It was found that extracting the average MFP from the fitting of the size-dependent effective thermal conductivity by a gray-medium model could well reflect the mode-dependent size effects. As shown in (12), the sum of the deviations with different channel thicknesses is defined as $\mathcal{L}(l_{\text{ave}})$. By employing Majumdar's model for cross-plane heat conduction [48] to calculate the effective thermal conductivity and minimize $\mathcal{L}(l_{\text{ave}})$, the average MFPs of different materials can be evaluated

$$\mathcal{L}(l_{\text{ave}}) = \left| \sum_t \frac{1}{3} \sum_j \int_0^{\omega_j} \hbar\omega \frac{\partial f_0}{\partial T} \text{DOS}_j(\omega) v_{g\omega_j} \frac{1}{1 + \frac{4}{3} \frac{l_t}{t}} d\omega - \frac{1}{1 + \frac{4}{3} \frac{l_{\text{ave}}}{t}} \right|^2 . \quad (12)$$

The extracted average MFPs of different materials are listed in Table II. With the average MFP, the average Knudsen numbers $Kn_t$ and $Kn_w$ can then be defined as (7) and (8). To clarify the $Kn_t$ dependence of the deviations between the predictions of the model and dispersion MC, a thermal resistance ratio $r_t = R_{\text{MC}}/R_{m,j,1}$ was introduced, where $R_{\text{MC}}$ and $R_{m,j,1}$ are the total thermal resistance predicted by dispersion MC and the model based on $l_{m,j,1}$, respectively. The thermal spreading effect and the ballistic effect with $w_g$ comparable with MFP are canceled in this ratio, thus, the ratio should reflect the cross-plane ballistic effect on the thermal resistance.
<table><caption>TABLE II
AVERAGE PHONON MFPs OF DIFFERENT SEMICONDUCTOR MATERIALS</caption><thead><tr><td>Material</td><td>Average MFP (nm)</td></tr></thead><tbody><tr><td>GaN</td><td>1612.3</td></tr><tr><td>AlN</td><td>3401.4</td></tr><tr><td>SiC</td><td>2506.9</td></tr><tr><td>β-Ga₂O₃</td><td>450.7</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Thermal resistance ratio $r_t$ as a function of $Kn_t$ of different materials.</figcaption></figure>

Fig. 5 shows the thermal resistance ratio varying with the average $Kn_t$ of different materials. It can be found that $r_t$ is always larger than 1 and increases with the increasing $Kn_t$, which more clearly shows that the model underestimates the cross-plane ballistic effect. It should be noted that in the model, the effective thermal conductivity is defined based on the total thermal resistance [24], that is, the average heat source temperature, which is more concerned in electronic devices. Whereas the way to involve the phonon dispersion with an integral over the phonon modes is originally used to consider the contribution of different phonon modes to the heat flux in the system with hot and cold phonon baths [23], [49]. Due to the nonlocal nature of the ballistic-diffusive heat transport [50], there is an inevitable deficiency in the gray model.

In addition, as shown in Fig. 5, the $r_t$ data of different materials collapse to a universal function by the average phonon MFP. The collapse essentially comes from the universal phonon MFP spectrum of these materials, which can be well depicted by Born-von Karman model [34]. Also, it can be seen that a logarithmic function can be used to depict the $r_t$ dependence of $Kn_t$

$$r_t = 0.15 \ln(Kn_t) + 1.35. \quad (13)$$

It should be noted that although the fit parameters are material independent, as discussed in [34], they are related to the phonon spectrum. At low temperatures, the phonon spectrum of different materials is markedly different from those at room

<figure><img src="image_6.png" /><figcaption>Fig. 6. Dimensionless total thermal resistance of GaN as a function of $w/t$ predicted by the model based on $l_{m,j,2}$ and MC simulations, with $w_g/w = 0.01$, $t = 0.2$, and $4 \ \mu m$, respectively.</figcaption></figure>

temperatures [51]. In this situation, the fitted parameters need to be adjusted.

### C. Ballistic Effect With $w_g$ Comparable With MFP

By introducing $r_t$ to the original modified MFPs, as shown in the following equation:

$$l_{m,j,2} = \frac{l_{0,j}}{(1 + A_w Kn_{w-\omega,j})(1 + \frac{2}{3} Kn_{t-\omega,j})r_t} \quad (14)$$

and calculating the effective thermal conductivity and the total thermal resistance using $l_{m,j,2}$, the differences between the model predictions and MC simulation results caused by the cross-plane ballistic effect can be compensated. Thus, the remaining deviations should be attributed to the ballistic effect with $w_g$ comparable with MFP.

Fig. 6 shows the dimensionless total thermal resistance of GaN predicted by dispersion MC and the model based on $l_{m,j,2}$. It can be seen that with $t = 4 \ \mu m$, the model predictions agree well with the MC simulation results. This indicates that the deviations with the larger $t$ mainly come from the model's underestimation of the cross-plane ballistic effect. In the case of $t = 0.2 \ \mu m$, when $w/t$ is large, the model's predictions coincide well with the MC simulation results, but the discrepancy increases significantly with decreasing $w/t$.

To clarify the $Kn_w$ dependence of this discrepancy, another thermal resistance ratio $r_{wg} = R_{MC}/R_{m,j,2}$ was introduced, where $R_{MC}$ and $R_{m,j,2}$ are the total thermal resistance predicted by dispersion MC and the model based on $l_{m,j,2}$, respectively. The thermal spreading effect and the cross-plane ballistic effect are canceled in this ratio, and thus, the ratio can reflect the influence of the ballistic effect with $w_g$ comparable with MFP on the thermal resistance.

Fig. 7 shows the thermal resistance ratio varying with the average $Kn_w$ of different materials. The $r_{wg}$ data of different materials also collapse to a universal function as $r_t$. It is found that $r_{wg}$ stays around unity with small $Kn_w$. This indicates that the fitted linear function $1 + A_w Kn_w$ in the gray model can well characterize the ballistic effect with $w_g$ comparable with MFP with small $Kn_w$. A logarithmic linear function can
<figure><img src="image_7.png" /><figcaption>Fig. 7. Thermal resistance ratio $r_{wg}$ as a function of $Kn_w$ of different materials, with $w/t = 2$ and 10.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8. $Kn_{w,0}$ as a function of $w_g/t$, with $w_g/w = 0.01, 0.02$, and 0.03.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Dimensionless total thermal resistance of different semiconductors as a function of $w/t$ predicted by the revised model and MC simulations, with $w_g/w = 0.01$ and (a) $t = 4$ µm, (b) $t = 0.2$ µm.</figcaption></figure>

be used to characterize the $Kn_w$ dependence of $r_{wg}$ after a threshold value $Kn_{w,0}$

$$r_{wg} = -0.17 \ln(Kn_w - Kn_{w,0}). \quad (15)$$

The slope of the function is nearly the same with different geometries, whereas the threshold value $Kn_{w,0}$ depends on the geometrical parameters since $A_w$ is a function of $w_g/w$ and $w/t$. From Fig. 8, it is found that $Kn_{w,0}$ is only a function of $w_g/t$. As $w_g/t$ increases, $Kn_{w,0}$ drops rapidly and approaches a plateau value. This implies that for most practical devices with a relatively large $w_g/t$, $Kn_{w,0}$ can be approximated with a constant.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Dimensionless temperature distributions of GaN along the heat source plane with $w_g/w = 0.01$ and $w/t = 40$.</figcaption></figure>

## D. Revised Thermal Resistance Model

We analyzed the deviations of the total thermal resistance predicted by dispersion MC and the semiempirical thermal resistance model. Two material-independent correction factors $r_t$ and $r_{wg}$ were introduced to the gray model

$$l_{m,j,r} = \frac{l_{0,j}}{(1 + A_w Kn_{w,0,j})(1 + \frac{2}{3} Kn_{t,0,j})r_t r_{wg}} \quad (16)$$

to compensate for its insufficiency to reflect the cross-plane ballistic effect and the ballistic effect with $w_g$ comparable with MFP, respectively. Although this work is focused on the thermal spreading process of WBG semiconductors, the same

analysis process can be extended to gray models for other systems [23] or to other materials.

Fig. 9 shows the dimensionless total thermal resistance predicted by dispersion MC and the revised model of different materials with $t = 4$ µm and $t = 0.2$ µm, respectively. The results of Fourier's law-based model are also given for comparison. As shown in Fig. 9, Fourier's law significantly underestimates the thermal resistance, which underscores the importance of taking the ballistic effects into account. Additionally, the separation of the curves of different materials
suggests that different effective thermal conductivities should be used to reflect the influence of phonon dispersion of different materials. The consistent agreement of the thermal resistance predicted by the revised model and dispersion MC verifies the effectiveness of the model.

**Fig. 10** further compares the dimensionless temperature distributions along the heat source plane. The dimensionless temperature, $θ$, is defined as $θ = (T - T_0)/(QR_{1d\_0})$. It more clearly shows that the use of bulk thermal conductivity significantly underestimates the junction temperature. The peak dimensionless temperature predicted by dispersion MC is nearly 2.7 times that from FEM with bulk thermal conductivity. Whereas FEM with the model predicted effective thermal conductivity can provide nearly the same dimensionless temperature distributions as the MC simulation.

## IV. CONCLUSION

We inspected the semiempirical thermal resistance model in a ballistic-diffusive regime using the phonon MC method with the phonon dispersion of several typical WBG semiconductor materials. This work shows that the gray model is insufficient to reflect the ballistic effects caused by the wide phonon MFP distributions. Two material-independent correction factors $r_t$ and $r_{wg}$ were introduced to the original model to improve its ability to reflect the cross-plane ballistic effect and the ballistic effect with $w_g$ comparable with MFP, respectively. The revised model can well address the thermal spreading process, ballistic effects, and the influence of phonon dispersion of different materials. Additionally, the model is easily coupled with widely used FEM-based thermal analysis. The results indicate that FEM with the model predicted effective thermal conductivity can give nearly the same junction temperature as the MC simulation. This work can provide a clearer understanding of the influence of phonon dispersion on the thermal transport process, and the model can be useful for the thermal management of WBG semiconductor devices. Additionally, the same analysis process and model construction methodology are expected to be applied in more complex and practical device situations, such as uneven internal heat source distributions [52], [53] and coupled with electro-thermal device simulations [54], [55].

## REFERENCES

[1] U. K. Mishra, P. Parikh, and Y.-F. Wu, “AlGaN/GaN HEMTs—An overview of device operation and applications,” *Proc. IEEE*, vol. 90, no. 6, pp. 1022–1031, Jun. 2002.

[2] F. Zeng et al., “A comprehensive review of recent progress on GaN high electron mobility transistors: Devices, fabrication and reliability,” *Electronics*, vol. 7, p. 377, Dec. 2018.

[3] S. Sriram, H. Hagleitner, D. Namishia, T. Alcorn, T. Smith, and B. Pulz, “High-gain SiC MESFETs using source-connected field plates,” *IEEE Electron Device Lett.*, vol. 30, no. 9, pp. 952–953, Sep. 2009.

[4] S. Zhu et al., “Improved MRD 4H-SiC MESFET with high power added efficiency,” *Micromachines*, vol. 10, no. 7, p. 479, Jul. 2019.

[5] M. Higashiwaki, K. Sasaki, A. Kuramata, T. Masui, and S. Yamakoshi, “Gallium oxide (Ga₂O₃) metal-semiconductor field-effect transistors on single-crystal β-Ga₂O₃ (010) substrates,” *Appl. Phys. Lett.*, vol. 100, no. 1, Jan. 2012, Art. no. 013504.

[6] M. H. Wong, K. Sasaki, A. Kuramata, S. Tamakoshi, and M. Higashiwaki, “Field-plated Ga₂O₃ MOSFETs with a breakdown voltage of over 750 V,” *IEEE Electron Device Lett.*, vol. 37, no. 2, pp. 212–215, Feb. 2016.

[7] C. Wang et al., “Progresses in state-of-the-art technologies of β-Ga₂O₃ devices,” *J. Phys. D, Appl. Phys.*, 2021.

[8] Y.-R. Wu and J. Singh, “Transient study of self-heating effects in AlGaN/GaN HFETs: Consequence of carrier velocities, temperature, and device performance,” *J. Appl. Phys.*, vol. 101, Jun. 2007, Art. no. 113712.

[9] Z. Cheng et al., “Significantly reduced thermal conductivity in β-(Al₀.₁Ga₀.₉)₂O₃/Ga₂O₃ superlattices,” *Appl. Phys. Lett.*, vol. 115, no. 9, 2019, Art. no. 092105.

[10] Y.-B. Liu, J.-Y. Yang, G.-M. Xin, L.-H. Liu, G. Csányi, and B.-Y. Cao, “Machine learning interatomic potential developed for molecular simulations on thermal properties of β-Ga₂O₃,” *J. Chem. Phys.*, vol. 153, no. 14, Oct. 2020, Art. no. 144501.

[11] P. Walterit et al., “Influence of AlGaN barrier thickness on electrical and device properties in Al₀.₁Ga₀.₈₆N/GaN high electron mobility transistor structures,” *J. Appl. Phys.*, vol. 112, no. 5, 2012, Art. no. 053718.

[12] K. R. Bagnall, “Device-level thermal analysis of GAN-based electronics,” Ph.D. dissertation, Dept. Mech. Eng., Massachusetts Inst. Technol., Cambridge, MA, USA, 2013.

[13] D.-S. Tang, G.-Z. Qin, M. Hu, and B.-Y. Cao, “Thermal transport properties of GaN with biaxial strain and electron-phonon coupling,” *J. Appl. Phys.*, vol. 127, no. 3, Jan. 2020, Art. no. 035102.

[14] D.-S. Tang and B.-Y. Cao, “Phonon thermal transport properties of GaN with symmetry-breaking and lattice deformation induced by the electric field,” *Int. J. Heat Mass Transf.*, vol. 179, Nov. 2021, Art. no. 121659.

[15] B. Chatterjee, K. Zeng, C. D. Nordquist, U. Singisetti, and S. Choi, “Device-level thermal management of gallium oxide field-effect transistors,” *IEEE Trans. Compon., Packag., Manuf. Technol.*, vol. 9, no. 12, pp. 2352–2365, Dec. 2019.

[16] B. Chatterjee et al., “Nanoscale electro-thermal interactions in AlGaN/GaN high electron mobility transistors,” *J. Appl. Phys.*, vol. 127, no. 4, Jan. 2020, Art. no. 044502.

[17] A. Sarua et al., “Thermal boundary resistance between GaN and substrate in AlGaN/GaN electronic devices,” *IEEE Trans. Electron Devices*, vol. 54, no. 12, pp. 3152–3158, Dec. 2007.

[18] M. Razavi, Y. S. Muzychka, and S. Kocabiyik, “Review of advances in thermal spreading resistance problems,” *J. Thermophys. Heat Transf.*, vol. 30, no. 4, pp. 863–879, Oct. 2016.

[19] Y. S. Muzychka, K. R. Bagnall, and E. N. Wang, “Thermal spreading resistance and heat source temperature in compound orthotropic systems with interfacial resistance,” *IEEE Trans. Compon., Packag., Manuf. Technol.*, vol. 3, no. 11, pp. 1826–1841, Nov. 2013.

[20] M. Garven and J. P. Calame, “Simulation and optimization of gate temperatures in GaN-on-SiC monolithic microwave integrated circuits,” *IEEE Trans. Compon. Packag. Technol.*, vol. 32, no. 1, pp. 63–72, 2009.

[21] A. A. Orouji and A. Aminbeidokhti, “A novel double-recessed 4H-SiC MESFET with partly undoped space region,” *Superlattices Microstructures*, vol. 50, no. 6, pp. 680–690, Dec. 2011.

[22] G. Chen, “Non-fourier phonon heat conduction at the microscale and nanoscale,” *Nature Rev. Phys.*, vol. 3, no. 8, pp. 555–569, 2021.

[23] Y.-C. Hua and B.-Y. Cao, “Ballistic-diffusive heat conduction in multiply-constrained nanostructures,” *Int. J. Therm. Sci.*, vol. 101, pp. 126–132, Mar. 2016.

[24] Y.-C. Hua, H.-L. Li, and B.-Y. Cao, “Thermal spreading resistance in ballistic-diffusive regime for GaN HEMTs,” *IEEE Trans. Electron Devices*, vol. 66, no. 8, pp. 3296–3301, Aug. 2019.

[25] H. Bao, J. Chen, X. Gu, and B.-Y. Cao, “A review of simulation methods in micro/nanoscale heat conduction,” *ES Energy Environ.*, vol. 1, pp. 16–55, Sep. 2018.

[26] T. Sadi and R. W. Kelsall, “Monte Carlo study of the electrothermal phenomenon in silicon-on-insulator and silicon-germanium-on-insulator metal-oxide field-effect transistors,” *J. Appl. Phys.*, vol. 107, no. 6, Mar. 2010, Art. no. 064506.

[27] Q. Hao, H. Zhao, and Y. Xiao, “A hybrid simulation technique for electrothermal studies of two-dimensional GaN-on-SiC high electron mobility transistors,” *J. Appl. Phys.*, vol. 121, no. 20, May 2017, Art. no. 204501.

[28] Q. Hao, H. Zhao, Y. Xiao, Q. Wang, and X. Wang, “Hybrid electrothermal simulation of a 3-D fin-shaped field-effect transistor based on GaN nanowires,” *IEEE Trans. Electron Devices*, vol. 65, no. 3, pp. 921–927, Mar. 2018.

[29] H.-L. Li, J. Shiomi, and B.-Y. Cao, “Ballistic-diffusive heat conduction in thin films by phonon Monte Carlo method: Gray medium approximation versus phonon dispersion,” *J. Heat Transf.*, vol. 142, no. 11, Nov. 2020, Art. no. 112502.
[30] A. Wang, M. J. Tadjer, and F. Calle, "Simulation of thermal management in AlGaN/GaN HEMTs with integrated diamond heat spreaders," *Semicond. Sci. Technol.*, vol. 28, no. 5, May 2013, Art. no. 055010.

[31] H. Zhang, Z. Guo, and Y. Lu, "Enhancement of hot spot cooling by capped diamond layer deposition for multifinger AlGaN/GaN HEMTs," *IEEE Trans. Electron Devices*, vol. 67, no. 1, pp. 47–52, Jan. 2020.

[32] C. Yuan et al., "Modeling and analysis for thermal management in gallium oxide field-effect transistors," *J. Appl. Phys.*, vol. 127, no. 15, Apr. 2020, Art. no. 154502.

[33] K. R. Bagnall, Y. S. Muzychka, and E. N. Wang, "Analytical solution for temperature rise in complex multilayer structures with discrete heat sources," *IEEE Trans. Compon., Packag., Manuf. Technol.*, vol. 4, no. 5, pp. 817–830, May 2014.

[34] J. P. Freedman, J. H. Leach, E. A. Preble, Z. Sitar, R. F. Davis, and J. A. Malen, "Universal phonon mean free path spectra in crystalline semiconductors at high temperature," *Sci. Rep.*, vol. 3, no. 1, pp. 1–6, Dec. 2013.

[35] Y. S. Muzychka, J. R. Culham, and M. M. Yovanovich, "Thermal spreading resistance of eccentric heat sources on rectangular flux channels," *J. Electron. Packag.*, vol. 125, no. 2, pp. 178–185, 2003.

[36] Y.-C. Hua and B.-Y. Cao, "Slip boundary conditions in ballistic-diffusive heat transport in nanostructures," *Nanosc. Microscale Thermophys. Eng.*, vol. 21, no. 3, pp. 159–176, Jul. 2017.

[37] D.-S. Tang, Y.-C. Hua, B.-D. Nie, and B.-Y. Cao, "Phonon wave propagation in ballistic-diffusive regime," *J. Appl. Phys.*, vol. 119, no. 12, 2016, Art. no. 124301.

[38] Y.-C. Hua and B.-Y. Cao, "Phonon ballistic-diffusive heat conduction in silicon nanofilms by Monte Carlo simulations," *Int. J. Heat Mass Transf.*, vol. 78, pp. 755–759, Nov. 2014.

[39] J.-P.-M. Péraud and N. G. Hadjiconstantinou, "Efficient simulation of multidimensional phonon transport using energy-based variance-reduced Monte Carlo formulations," *Phys. Rev. B, Condens. Matter*, vol. 84, no. 20, Nov. 2011, Art. no. 205331.

[40] J. D. Chung, A. J. H. McGaughey, and M. Kaviani, "Role of phonon dispersion in lattice thermal conductivity modeling," *J. Heat Transf.*, vol. 126, no. 3, pp. 376–380, Jun. 2004.

[41] G. Chen, *Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons*. London, U.K.: Oxford Univ. Press, 2005.

[42] Z. Guo et al., "Anisotropic thermal conductivity in single crystal β-gallium oxide," *Appl. Phys. Lett.*, vol. 106, no. 11, Mar. 2015, Art. no. 111909.

[43] Z. Yan and S. Kumar, "Phonon mode contributions to thermal conductivity of pristine and defective β-Ga₂O₃," *Phys. Chem. Chem. Phys.*, vol. 20, no. 46, pp. 29236–29242, 2018.

[44] P. Jiang, X. Qian, X. Li, and R. Yang, "Three-dimensional anisotropic thermal conductivity tensor of single crystalline β-Ga₂O₃," *Appl. Phys. Lett.*, vol. 113, no. 23, Dec. 2018, Art. no. 232105.

[45] B. A. Danilchenko, I. A. Obukhov, T. Paszkiewicz, S. Wolski, and A. Jeżowski, "On the upper limit of thermal conductivity GaN crystals," *Solid State Commun.*, vol. 144, nos. 3–4, pp. 114–117, Oct. 2007.

[46] G. A. Slack, R. A. Tanzilli, R. Pohl, and J. Vandersande, "The intrinsic thermal conductivity of ALN," *J. Phys. Chem. Solids*, vol. 48, no. 7, pp. 641–647, 1987.

[47] Q. Zheng, C. Li, A. Rai, J. H. Leach, D. A. Broido, and D. G. Cahill, "Thermal conductivity of GaN, 71GaN, and SiC from 150 K to 850 K," *Phys. Rev. Mater.*, vol. 3, no. 1, Jan. 2019, Art. no. 014601.

[48] A. Majumdar, "Microscale heat conduction in dielectric thin films," *J. Heat Transf.*, vol. 115, no. 1, pp. 7–16, Feb. 1993.

[49] N. Mingo, "Calculation of Si nanowire thermal conductivity using complete phonon dispersion relations," *Phys. Rev. B, Condens. Matter*, vol. 68, no. 11, Sep. 2003, Art. no. 113308.

[50] H.-L. Li and B.-Y. Cao, "Radial ballistic-diffusive heat conduction in nanoscale," *Nanosc. Microscale Thermophys. Eng.*, vol. 23, no. 1, pp. 10–24, Jan. 2019.

[51] K. T. Regner, D. P. Sellan, Z. Su, C. H. Amon, A. J. H. McGaughey, and J. A. Malen, "Broadband phonon mean free path contributions to thermal conductivity measured using frequency domain thermoreflectance," *Nature Commun.*, vol. 4, no. 1, pp. 1–7, Jun. 2013.

[52] X. Chen, S. Boumaiza, and L. Wei, "Self-heating and equivalent channel temperature in short gate length GaN HEMTs," *IEEE Trans. Electron Devices*, vol. 66, no. 9, pp. 3748–3755, Sep. 2019.

[53] Y.-C. Hua and B.-Y. Cao, "The effective thermal conductivity of ballistic-diffusive heat conduction in nanostructures with internal heat source," *Int. J. Heat Mass Trans.*, vol. 92, pp. 995–1003, Jan. 2016.

[54] X. Chen, S. Boumaiza, and L. Wei, "Modeling bias dependence of self-heating in GaN HEMTs using two heat sources," *IEEE Trans. Electron Devices*, vol. 67, no. 8, pp. 3082–3087, Aug. 2020.

[55] S. Choi et al., "A perspective on the electro-thermal co-design of ultra-wide bandgap lateral devices," *Appl. Phys. Lett.*, vol. 119, no. 17, Oct. 2021, Art. no. 170501.