Modeling effective thermal conductivity enhanced by surface waves using the Boltzmann transport equation

The thermal management of semiconductors at the device level has become a crucial issue owing to the high integration density and miniaturization of microelectronic systems. Because surface phonon polaritons (SPhPs) exhibit long propagation lengths, they are expected to contribute significantly to the heat dissipation in microelectronic systems. This study aims to numerically estimate the heat transfer due to SPhPs in a thin SiO2 film. The one-dimensional Boltzmann transport equation (BTE) is solved using the estimated propagation length based on the SPhP dispersion curves. The temperature profiles and heat fluxes are predicted and demonstrate the size effect of the film on the effective in-plane thermal conductivity of the SiO2 film. The results indicate that the temperature distribution was constant regardless of the film length and thickness because the propagation length was much longer than the film length. In addition, the heat flux increased with decreasing film thickness owing to the depth-averaged energy transfer. The effective thermal conductivities predicted using the BTE differed by ~ 16.5% from the values obtained from the analytical expression. The numerical results of this study can provide valuable data when studying the thermal behavior of SPhPs.

Surface phonon polaritons (SPhPs) are the collective motion of optical phonons coupled with photons (or electromagnetic waves) in polar dielectrics. Physically, SPhPs can have long propagation lengths of up to several millimeters or even several meters, which is several orders of magnitude greater than their phonon mean free path, along the surface of polar dielectrics 1-3 . Therefore, SPhPs can ballistically transport thermal energy. Many studies have investigated SPhPs as additional thermal energy carriers owing to their long propagation lengths [4][5][6][7] . Two decades ago, Chen et al. 8 first predicted the effective in-plane thermal conductivity of SPhPs (κ SPhPs ) in a freestanding SiO 2 film that is isotropic and homogeneous materials. Their results indicate that κ SPhPs increases dramatically when the film thickness is less than 100 nm. Recently, Tranchant et al. 9 and Wu et al. 10 experimentally demonstrated the variations in κ SPhPs with the thickness, length, and temperature of a thin film. Owing to their enhanced effective in-plane thermal conductivity, SPhPs are potentially useful for heat dissipation purposes in next-generation microelectronics systems. Currently, the thermal management in microelectronics is mainly conducted at the package level. Nevertheless, device-level thermal management has garnered substantial attention because of their high integration density and miniaturization of microelectronic systems 11,12 . Therefore, predicting the heat transfer capabilities of SPhPs is important; however, their temperature distribution and heat flux have rarely been analyzed. In this paper, the Boltzmann transport equation (BTE) is used instead of Fourier's law to predict the temperature profile and heat flux at the nanoscopic scale. It is well known that Fourier's law cannot predict heat conduction at the submicron scale, where the characteristic length scale is comparable to, or smaller than, the mean free path of the energy carrier [13][14][15][16] . Conversely, BTE can describe the non-diffusive (i.e., ballistic) heat conduction at such a small scale and can consider spatial behavior in anisotropic and inhomogeneous materials or film defects. For example, Majumdar 17 and Joshi and Majumdar 18 predicted the temperature distribution within a nanoscale diamond substrate. Specifically, they derived the equation of phonon radiative transfer in terms of the phonon intensity from the BTE using the relaxation time approximation. Chen et al. [19][20][21] developed a ballistic-diffusive equation (BDE) based on the BTE to describe the diffusive and ballistic energy transport contributions of the energy carriers and calculated the temperature and heat flux in the imaginary material. These studies show that the BTE can describe nanoscale heat transfer well.
This study aims to quantitatively analyze the heat transfer due to SPhPs in a nanoscale thin film using the BTE and lay the foundation for heat dissipation design considering SPhPs. Specifically, the one-dimensional BTE is established and directly solved to predict the temperature field, heat flux, and effective in-plane thermal conductivity of SPhPs in a freestanding SiO 2 film. The effects of the length and thickness of the SiO 2 film on the effective in-plane thermal conductivity of the SPhPs are further examined using the established BTE. The numerical solution to the BTE for the SPhPs in a freestanding SiO 2 film is compared with the theoretical value obtained by Tranchant et al. 9 Mathematical modeling Boltzmann transport equation for surface phonon polaritons. The one-dimensional BTE is expressed with the relaxation time approximation 22 as where f is the distribution function of the SPhPs, v g is the group velocity of the SPhPs, f 0 is the equilibrium distribution function (Bose-Einstein distribution function), and τ is the relaxation time. Assuming that the film is sufficiently thin such that the heat transfer in the thickness direction can be neglected, the one-dimensional steady-state BTE along the in-plane direction (x-direction in Fig. 1) can be rewritten using the intensity notation as where Λ e is the effective propagation length of the SPhP obtained using Matthiessen's rule; that is, where Λ SPhP is the propagation length of the SPhP and L is the film length 1 . In addition, θ is the polar angle along the x-axis, and I ω and I 0 ω are the directional SPhP intensity for a particular frequency per www.nature.com/scientificreports/ unit length (W s m −1 rad −1 ) and the equilibrium SPhP intensity, respectively, and are defined as (Supplementary information S1) where ħ is the Planck constant divided by 2π, ω is the SPhP frequency, and D 2D (ω) is the SPhP density of states per unit area. Once the SPhP intensity is determined, the temperature distribution and heat flux (q SPhPs ) can be obtained by integrating over the SPhP frequency range as where ω H and ω L are the highest and lowest frequencies of the SPhPs, respectively, and d is the film thickness. In Eq. (6), the heat flux is inversely proportional to the film thickness. The heat flux of a surface polariton includes the assumption that the amount of in-plane energy transferred by the SPhPs is constant regardless of the position according to film thickness direction and is a depth-averaged value. Because SPhPs can propagate in any direction on the surface of a thin film, the numerical method used to solve the BTE should include the discretization of space and polar angle. In the present study, the finite element method was used with the discrete ordinates method (DOM) to discretize the spatial and angular terms of the BTE. The DOM numerically integrates angles by assigning a weight function w to a specific angle 23 . The Gaussian-Legendre quadrature was adopted to integrate in polar angle. The boundary conditions of the BTE were given by the constant equilibrium intensity of the SPhPs at a fixed temperature at the boundary. Because the SPhPs propagate forward (positive x-direction) when cosθ > 0 and backward (negative x-direction) when cosθ < 0, the boundary conditions for the two cases can be expressed as 18 where T l and T r are the temperatures at the boundaries of the analysis domain, as shown in Fig. 1.
Effective in-plane thermal conductivity. The effective in-plane thermal conductivity due to the SPhPs in the layered thin film is given as follows 8,9 : where β R is the real part of the SPhP in-plane wave vector β. Equation (8) expresses the effective in-plane thermal conductivity at an interface of infinite length. The effective propagation length is used in consideration of the boundary scattering caused by the propagation length of the SPhPs, which is significantly longer than the film length 1,9 .
The optical phonons are thermally excited in the frequency range of 7.6-258 Trad/s for SiO 2 9 . In this frequency range, the SPhPs can be classified into three modes: Zenneck, SPhP, and transverse magnetic (TM)-guided modes 24 . However, the TM-guided mode was excluded because it does not satisfy the existence condition of SPhPs on a SiO 2 film 9 . Therefore, we considered both the SPhP and Zenneck modes in this study.
Dispersion relation of surface phonon polaritons. The dispersion relation estimates the physical properties of energy carriers such as propagation length and group velocity as a relationship between frequency and wave vector. These properties were required to solve the BTE and obtain theoretical solutions. www.nature.com/scientificreports/ The dispersion relation for an amorphous SiO 2 thin film, assumed to be a nonmagnetic material surrounded by air as shown in Fig. 2, can be derived from Maxwell's equations for the three layers of the thin film 25 : where p j is the transverse wave vector of the SPhPs given by p 2 j = β 2 − ε j k 2 0 , ε j is the relative permittivity of the material (j = 0, 1, 2), and k 0 is the wave vector defined as the frequency divided by the speed of light in a vacuum.
For lossy dielectric materials, the permittivity can have a complex form and can be expressed as the relationship between the refractive index n and the extinction coefficient k of the material; that is, ε = (n + i•k) 2 , where i is an imaginary number. The permittivity of amorphous SiO 2 at room temperature (300 K) was obtained from the experimental data, using its refractive index and extinction coefficient 26 .

Results and discussion
To consider the propagating nature of the SPhPs along the surface of the SiO 2 film, the wave vector of the SPhPs was set to a complex number (i.e., β = β R + i•β I ) and the frequency was set to a real number. From the dispersion relation, the propagation length of the SPhPs was defined as Λ SPhP = 1/(2•β I ) 27 .
The dispersion relation consists of the frequency and the real part of the wave vector, as shown in Fig. 3a. Because the trend of the dispersion relation describing the SPhPs in SiO 2 films of different thicknesses aligns with the light line, it can be deduced that the group velocity (v g = dω/dβ R ) of the SPhPs is comparable to the speed of light. Thus, the SPhP density of states per unit area can be expressed as   www.nature.com/scientificreports/ Figure 3b shows that the propagation length of the SPhPs on the surface of a freestanding SiO 2 film increases with decreasing film thickness. In addition, when the film thickness was 25 nm, the SPhPs propagated for up to several meters, particularly at low frequencies. Therefore, the long propagation length of SPhPs is expected to result in high thermal conductivity for thinner films.
In this study, the same SiO 2 film dimensions used in the report by Tranchant et al. 9 were applied as a preliminary study to compare the experimental and theoretical results of the effective in-plane thermal conductivity of SPhPs. Thus, the film lengths were 120 μm, 195 μm, 275 μm, and 350 μm, and the film thicknesses were 25 nm, 50 nm, 75 nm, and 200 nm. The BTE was solved for 16 cases with each length and thickness of the film. For all cases, the temperature difference to generate heat transfer was set to ∆T = T l -T r = 1 K, where T r was set to 300 K. The BTE was solved using COMSOL Multiphysics 5.6, a commercial package that was employed to compute the BTE in several studies 28,29 .
Grid independence tests were performed for both the angular and spatial discretization for a film length of 120 μm and film thickness of 25 nm by monitoring the heat flux and temperature difference between both ends of the film. The number of Gaussian points and elements ranged between 4-32 and 20-60, respectively. The temperature difference and heat flux gradually converged as the number of Gaussian points increased. The deviations in the temperature difference and heat flux between the 16th and 32nd Gaussian points were 0.077% and 0.053%, respectively. Hence, in this study, the number of Gaussian points was set to 16. Conversely, even if the number of elements increased, the computational time and temperature difference remained; therefore, it was decided that 60 elements would be used. Figure 4 shows the temperature distribution of the freestanding SiO 2 film for various film lengths and thicknesses. Temperature jumps caused by the ballistic energy transport of the SPhPs were observed at the film boundary for all cases. The temperature distributions were almost identical regardless of the film length and thickness; however, the temperature gradient (∆T/L) decreased with an increase in the film length, where ∆T is the temperature difference. This uniformity is due to the effective propagation length used in the BTE. Since the propagation length of the SPhPs on the SiO 2 film was significantly greater than the film length, as shown in Fig. 3b, the effective propagation length using Mathiessen's rule became equal to the film length for most frequencies (Λ e~L ). Therefore, the SPhPs exhibited the same behavior as that of boundary scattering for all film-length ranges. For this reason, we expect that there would be temperature differences when the film length is greater than the propagation length for most frequency ranges. Figure 5 shows the effect of the SiO 2 film thickness on the heat flux of the SPhPs for various film lengths. The heat flux was also nearly the same regardless of the film length, similar to the temperature distribution. However, the heat flux increased with a decrease in the film thickness and reached a maximum value at d = 25 nm. The effective propagation length was almost equal to the film length in this study. Therefore, the divergences in the SPhPs intensities calculated by the BTE based on the film length were not large. The heat flux was more affected by the thickness of the film, as described by Eq. (6). The heat flux may differ based on the film length when the effective propagation lengths of the SPhPs are not equal to the film length for most frequency ranges.
The effective in-plane thermal conductivity is defined as κ BTE = q SPhPs •L/∆T, based on Fourier's law. We predicted the effective in-plane thermal conductivity using the temperature distribution and heat flux calculated from the numerical solution to the BTE, as shown in Fig. 6 (solid line). The effective in-plane thermal conductivity increased as the film thickness decreased and the length increased. This is because the gradient of the temperature distribution decreases with increasing film length and the heat flux increases with decreasing Zenneck modes contribute more to the enhancement of the effective thermal conductivity among surface waves as reported in a previous study 9 . Since the same frequency range was used in the numerical method, the effective thermal conductivity was also enhanced by the Zenneck modes. The enhancement in the thermal conductivity clearly indicates that the surface waves at both the SiO 2 -air interfaces can be useful for heat dissipation purposes in microelectronics. Figure 6 also compares the calculated effective in-plane thermal conductivity with the theoretical values (see Eq. (8)). The effective in-plane thermal conductivity calculated by the BTE exhibits the same tendency as that of the values obtained from theory, but its value is always less by ~16.5% for all cases. It is worth noting that we directly solved the original form of the BTE in Eq. (2). In contrast, Tranchant et al. 9 considered that only a diffusive part of the SPhPs would contribute to the thermal conductivity; that is, they calculated the effective in-plane thermal conductivity based on the diffusion approximation, ∂f/∂x ≈ ∂f 0 /∂x.  www.nature.com/scientificreports/

Conclusions
In this study, the BTE for SPhPs was numerically solved to investigate the ballistic energy transport due to the SPhPs on the surface of an amorphous SiO 2 thin film. In particular, the dispersion relation of SPhPs on an amorphous SiO 2 film was used to obtain the propagation length and group velocity of the SPhPs for different film thicknesses and lengths. The film temperature distribution and heat flux were estimated based on the SPhP intensity. The conclusions of this study are as follows: 1. From the dispersion relation, it was found that the SPhPs on the amorphous SiO 2 film surrounded by air had long propagation lengths of several meters. It was confirmed that the SPhPs with long propagation lengths were the energy carriers that were conducive to ballistic energy transport; 2. The temperature distribution was almost identical regardless of the film length and thickness because the effective propagation length was equal to the film length for most frequency ranges and exhibited temperature jumps owing to the ballistic energy transport due to the SPhPs. However, the heat flux increased with the decrease in the film thickness owing to the depth-averaged energy transfer; 3. The effective in-plane thermal conductivity of the SPhPs increased with increasing film length and decreasing film thickness. In addition, the effective thermal conductivity calculated by the BTE was ~ 16.5% lower than the theoretical results. The current numerical approach can provide useful information regarding the thermal characteristics of SPhPs.
For further studies, the BTE should be solved for conditions under which the film length is greater than the propagation length of the SPhPs. In addition, the numerical solution to the BTE can be readily extended to 2-D simulations, and SPhP-enhanced thermal conductivity in a more complicated geometry can thus be considered by the proposed method. The prediction of the temperature distribution and effective in-plane thermal conductivity due to the SPhPs is expected to be useful in designing various products, such as semiconductors, in microelectronics.