Phonon Confinement and Transport in Ultrathin Films

TL;DR

This study investigates phonon-mediated thermal transport in ultrathin films using advanced lattice dynamics and molecular dynamics simulations, revealing size-dependent behaviors and confinement effects that challenge traditional models.

Contribution

It introduces a two-dimensional lattice dynamics approach that accurately captures phonon confinement effects in ultrathin films, improving upon traditional three-dimensional models.

Findings

Thermal conductivity decreases with film thickness down to 2-3 nm.

Confinement effects emerge below 5-6 nm thickness.

Bulk phonon properties with boundary scattering models fail at low thicknesses.

Abstract

Thermal transport by phonons in films with thicknesses of less than 10 nm is investigated in a soft system (Lennard-Jones argon) and a stiff system (Tersoff silicon) using two-dimensional lattice dynamics calculations and the Boltzmann transport equation. This approach uses a unit cell that spans the film thickness, which removes approximations related to the finite cross-plane dimension required in typical three-dimensional-based approaches. Molecular dynamics simulations are performed to obtain finite-temperature structures for the lattice dynamics calculations and to predict thermal conductivity benchmarks. Thermal conductivity decreases with decreasing film thickness until the thickness reaches four unit cells (2.1 nm) for argon and three unit cells (1.6 nm) for silicon. With a further decrease in film thickness, thermal conductivity plateaus in argon while it increases in silicon. This unexpected behavior, which we identify as a signature of phonon confinement, is a result of an increased contribution from low-frequency phonons, whose density of states increases as the film thickness decreases. Phonon mode-level analysis suggests that confinement effects emerge below thicknesses of ten unit cells (5.3 nm) for argon and six unit cells (3.2 nm) for silicon. Thermal conductivity predictions based on the bulk phonon properties combined with a boundary scattering model do not capture the low thickness behavior. To match the two-dimensional lattice dynamics and molecular dynamics predictions for larger thicknesses, the three-dimensional lattice dynamics calculations require a finite specularity parameter that in some cases approaches unity. These findings point to the challenges associated with interpreting experimental thermal conductivity measurements of ultrathin silicon films, where surface roughness and a native oxide layer impact phonon transport.

Tables

Table 1. Table 1: Lattice constants, linear RMS displacements along one of the Cartesian directions, and thermal conductivities ( k MD − GK subscript 𝑘 MD GK k_{\rm MD-GK} and k 3 ​ D − BTE subscript 𝑘 3 D BTE k_{\rm 3D-BTE} ) for bulk argon at temperatures of 5, 10, and 20 K and bulk silicon at a temperature of 200 K. The 3D-BTE results are converged to within 1% based on the wave vector resolution.

System	Argon			Silicon
Temperature	5 K	10 K	20 K	200 K
Lattice Constant (nm)	0.5269	0.5279	0.5301	0.5440
Linear RMS Displacement (nm)	0.0048	0.0069	0.0101	0.0054
$k_{\mathrm{MD}-\mathrm{GK}}$ (W/m-K)	$8.9\pm 0.2$	$3.9\pm 0.2$	$1.5\pm 0.1$	$414\pm 7$
$k_{3D-\mathrm{BTE}}$ (W/m-K)	9.0	4.1	1.7	464

Equations

$$\varphi (r_{ij})=4{ϵ}_{\mathrm{LJ}}\left[{\left(\frac{{\sigma }_{\mathrm{LJ}}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{\mathrm{LJ}}}{r_{ij}}\right)}^6\right].$$

$$k_{\alpha }=\sum _{\mathbf{\kappa },\nu }{c}_{\mathrm{ph}}\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)v_{\mathrm{g},\alpha }^2\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right){\tau }_{\alpha }\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right),$$

$${\tau }_{\alpha }\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)=F\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right){\tau }_{\mathrm{pp},\alpha }\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right),$$

$$F\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)=1-\frac{1-p\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)}{\delta \text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)}\frac{1-\exp [-\delta \text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)]}{1-p\text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)\exp [-\delta \text{mspace}-4.0mu\left(\text{mspace}-8.0mu\begin{array}{cc}& \mathbf{\kappa }\\ & \nu \end{array}\text{mspace}-3.0mu\right)]}.$$

$${\kappa }_z=\frac{2\pi l_z}{a{N}_z},$$

Full text

††preprint:

I Introduction

Advances in nano-fabrication have pushed the minimum feature size in electronic, optoelectronic, and energy conversion devices to below 10 nm.Heremans et al. (2013); Ieong et al. (2004); Wang et al. (2012) Thin films are prototypical nanostructures that are used in thermoelectric energy conversion, Venkatasubramanian et al. (2001) solar energy conversion, Wang et al. (2014) semiconductor laser,Lu et al. (2012) and solid state lightingYing et al. (2014) devices. The thermal transport properties of a film can vary dramatically compared to its corresponding bulk material.Cahill et al. (2003, 2014) In an electrically-insulating or semiconducting material, the dominate heat carriers are phonons, which scatter more often with a film’s boundaries as its thickness is reduced. The resulting thermal conductivity reduction is problematic when heat dissipation is important (e.g., in a laser), but desirable in thermoelectric energy conversion.

Extensive experimental measurements have been performed on films where phonons dominate thermal transport, demonstrating the significant in-plane thermal conductivity reduction as the film thickness is reduced. A collection of data for silicon films at a temperature of 300 K with thicknesses between 9 nm and 1.5 $\mu$m is plotted in Fig. 1. The measurements were performed using $\rm 3-\omega$,Asheghi et al. (1998); Ju and Goodson (1999); Ju (2005); Liu and Asheghi (2006); Hao et al. (2006) thermoreflectance (TR),Aubain and Bandaru (2010, 2011) suspended island,Tang et al. (2010); Yu et al. (2010) Raman thermometry,Chávez-Ángel et al. (2014); Neogi et al. (2015) and transient thermal grating (TTG)Cuffe et al. (2015) techniques. For the $\rm 3-\omega$ and TR measurements, the films were embedded between a metal capping layer and a substrate. For the suspended island, Raman thermometry, and TTG measurements, the films were suspended. While the thickness-dependent trend is clear, the thermal conductivities at any given thickness are scattered, with a spread of up to 48 W/m-K at a thickness of 150 nm, where the average value is 72 W/m-K. The same thermal conductivity trend, i.e., a reduction with decreasing film thickness, has been measured for GaN,Beechem et al. (2016) black phosphorus,Luo et al. (2015) and h-BN.Jo et al. (2013)

The origin of the thermal conductivity reduction in a film can be probed theoretically. Phonon-boundary scattering reduces the bulk mean free paths (MFPs) and thus reduces thermal conductivity. The film MFP for a given phonon mode can be obtained by combining the intrinsic bulk phonon-phonon scattering MFP with a boundary scattering model using the Matthiessen ruleMcGaughey et al. (2011); Jain et al. (2013); Wang and Huang (2014) or the Fuchs-Sondheimer (FS) suppression function.Cuffe et al. (2015) Phonon-boundary scattering can also be simulated using numerical techniques,McGaughey and Jain (2012); Yang and Chen (2004) which are useful for arbitrary geometries such as nanoporous films and nanocomposites.Jain et al. (2013); Tang et al. (2013); Fu et al. (2014, 2017) All of these methods use bulk phonon properties as the starting point. The phonon transport is thus three-dimensional (3D) and is based on the bulk unit cell shown in Fig. 2, which has a corresponding 3D Brillouin zone. While the focus here is on in-plane thermal transport, similar tools can be applied to investigate cross-plane thermal transport.Sellan et al. (2010); Vermeersch et al. (2016)

The results from a series of in-plane thermal conductivity predictions based on the 3D treatment for silicon are plotted in Fig. 1.McGaughey et al. (2011); Jain et al. (2013); Wang and Huang (2014); Fu et al. (2017); Cuffe et al. (2015) The bulk phonon properties were obtained from lattice dynamics calculations and a solution of the Boltzmann transport equation (BTE), which we will refer to as 3D-BTE, under the relaxation time approximation.McGaughey et al. (2011); Jain et al. (2013); Wang and Huang (2014); Fu et al. (2017); Cuffe et al. (2015) In the lattice dynamics calculations, the atomic interactions were modeled using an empirical modelMcGaughey et al. (2011); Cuffe et al. (2015); Holland (1963) or density functional theory calculations.Jain et al. (2013); Wang and Huang (2014); Fu et al. (2017); Cuffe et al. (2015); Esfarjani et al. (2011) All calculations were performed on suspended films with an assumption of fully-diffuse phonon-boundary scattering. As with the experimental measurements, the thermal conductivity predictions at any given film thickness are scattered, with a spread of up to 27 W/m-K at a thickness of 1 $\rm\mu$m, where the average value is 107 W/m-K.

The agreement between the measurements and predictions plotted in Fig. 1 worsens as the film thickness gets smaller. Differences may originate from: (i) Some of the experimental films being embedded, while all the predictions are for suspended films. For an embedded film, there is the potential for film phonons to transmit into the surrounding materials and vice-versa, while for a suspended film, all phonon-boundary interactions lead a specular or diffuse reflection. (ii) The presence of a native oxide layer on the experimental films.Neogi et al. (2015) While thermal transport in the oxide layer is not considered explicitly in the predictions, a predetermined parameter (i.e., the specularity) can be used to model the surface roughness. (iii) The validity of the theoretical framework. Our objective here is to examine (iii). Specifically, we will probe the assumption of using phonon modes based on the 3D bulk description.

The rigorous description of the phonons modes in a film is based on a unit cell that spans the film thickness and has a corresponding two-dimensional (2D) Brillouin zone, as shown in Fig. 2. Phonons only advect in the in-plane directions and the concepts of cross-plane transport and boundary scattering are not meaningful. For thick enough films, as shown in Fig. 1, the phonon transport can be adequately described using the 3D Brillouin zone and a boundary scattering model. But for thin enough films, which we will refer to as ultrathin films, the 3D treatment breaks down. This phenomenon is known as phonon confinement.Turney et al. (2010); Balandin and Wang (1998) The identification of the phonon confinement regime is a central objective of this work.

Kim et al.Kim et al. (2019) recently proposed an algorithm to map between 2D and 3D phonon modes, which they applied to Lennard-Jones (LJ) argon and multi-layer graphene films. A near one-to-one correspondence was discovered for thicknesses greater than 5 nm for LJ argon and for all thicknesses of graphene. Their study only considered harmonic phonon properties, however, and not the impact of confinement on phonon scattering and thermal conductivity. 2D harmonic lattice dynamics calculations have also been applied to predict the phonon properties of silicon films with surface nanopillars.Davis and Hussein (2014); Honarvar and Hussein (2018)

Beyond harmonic lattice dynamics, the 2D treatment is commonly applied to model phonon transport and scattering in 2D materials such as graphene,Lindsay et al. (2010a) multi-layer graphene,Lindsay et al. (2011) Mo$\rm S_{2}$,Li et al. (2013) and black phosphoreneJain and McGaughey (2015) using anharmonic lattice dynamics calculations and the BTE. While the 2D treatment is rigorous, the computational cost of the BTE-based method, which scales with the number of atoms in the unit cell to the fourth power, makes calculations on films challenging. To our knowledge, the 2D treatment at the anharmonic level has not been applied to films built from common semiconductors such as silicon or GaN.

Instead, previous researchers investigated confinement using the bulk 3D phonon modes. In the phonon depletion approach, only bulk phonon modes with cross-plane wavelengths that fit inside the film are considered.Turney et al. (2010); Wang and Huang (2014); Lindsay et al. (2010b) I.e., the finite sampling of the 3D Brillouin zone in the cross-plane direction is set by the film thickness. Phonon-boundary scattering must still be included. Turney et al.Turney et al. (2010) used the phonon depletion approach in 3D-BTE and the Green-Kubo (GK) method in MD simulations to predict the thermal conductivity of Lennard-Jones (LJ) argon films with thickness as small as 1 nm and Stillinger-Weber silicon films with thickness as small as 17.4 nm. Wang and HuangWang and Huang (2014) also investigated in-plane thermal transport in silicon films using 3D-BTE and phonon depletion with input from density functional theory calculations.

Herein, we perform 2D lattice dynamics calculations combined with an iterative solution to the BTE (2D-BTE) to study phonon transport in suspended films of LJ argon (a soft material) and Tersoff silicon (a stiff material) with thickness as small as two conventional unit cells (i.e., four atomic layers for argon and eight atomic layers for silicon). The 2D-BTE approach uses a unit cell that spans the film thickness (Fig. 2), thus directly treating the phonon transport as 2D without any phonon-boundary scattering. In doing so, we remove all approximations associated with the finite cross-plane dimension used in previous 3D-BTE studies. MD simulations are performed to obtain the finite temperature atomic structures and to predict thermal conductivity benchmarks. We also perform 3D-BTE calculations with the Fuchs-Sondheimer suppression function (3D-BTE-FS) and phonon-depletion to assess their validity. This multi-method approach allows us to unravel the similarities and differences between phonon transport in corresponding 2D film and 3D bulk systems. We use empirical potentials to describe the atomic interactions due to the large computational costs associated with the large unit cells. Our findings will be relevant to other film systems.

The rest of the paper is organized as follows. The potential parameters, brief descriptions of the MD and BTE methods, and computational details are provided in Sec. II. Results and discussion for the LJ argon and silicon films are presented in Secs. III and IV. Included are descriptions of the finite temperature atomic structures, anharmonicity analysis performed using root-mean-squared (RMS) displacements, and phonon properties and thermal conductivities predicted from the 2D-BTE, 3D-BTE, and MD treatments. Of note are predictions of unusual thermal conductivity trends for very thin films, underestimations of thermal conductivity from the 3D-BTE-FS description, and incorrect trends predicted by the depletion approach. The results are interpreted using the density of states (DOS), mode lifetimes, and thermal conductivity accumulation functions. The findings are summarized in Sec. V.

II Methodology

II.1 Interatomic Potentials and Structures

Solid argon is a face-centered cubic crystal with a one-atom basis. The potential energy $\phi$ between two atoms $i$ and $j$ is modeled using the two-body LJ potentialAshcroft and Mermin (1976)

[TABLE]

Here, $r_{ij}$ is the distance between the two atoms, $\epsilon_{\rm LJ}$ is the energy scale ($1.67\times 10^{-21}$ J), and $\sigma_{\rm LJ}$ is the length scale ($3.4\times 10^{-10}$ m). The atomic mass is $6.63\times 10^{-26}$ kg. The potential energy is shifted at a cut-off radius of $3\sigma_{\rm LJ}$. We found that using the typical cutoff radius of $2.5\sigma_{\rm LJ}$ led to imaginary modes in the film phonon dispersion curves that disappeared when the cutoff radius was increased to $3\sigma_{\rm LJ}$. Temperatures of 5, 10, and 20 K are considered. These low temperatures (the melting temperature of LJ argon is around 90 K) allow for comparisons with the BTE predictions, which include up to three-phonon scattering events. The significant anharmonicity in LJ argon at higher temperatures leads to non-negligible four-phonon scattering, which is not considered here.Turney et al. (2009); Feng and Ruan (2016)

Silicon has a diamond crystal structure and is modeled using the three-body Tersoff potential.Tersoff (1989) All calculations are performed at a temperature of 200 K to limit the effects of four-phonon scattering.Feng and Ruan (2016) We do not consider any reconstruction of the film surfaces. All structures are stable in the MD simulations and exhibit only real frequencies in their phonon dispersions.

The films are initialized as supercells built from $N_{x}\times N_{y}\times N_{z}$ bulk conventional unit cells (i.e., a simple cubic lattice with a four-atom basis for argon and an eight-atom basis for silicon). An in-plane size of $N_{x}=N_{y}=8$ is used, while $N_{z}$ is the number of unit cells in the cross-plane direction. Periodic boundary conditions are applied in all three directions. A vacuum region of three conventional units cells, which is larger than the range of the atomic interactions, is placed between the top and bottom surfaces of the film to ensure that the film is isolated and suspended.

Computational costs limit the film thicknesses that can be explored, even with the use of empirical potentials. For an efficient MD simulation, the computational cost scales linearly with the number of atoms in the computational cell. For the 2D-BTE and 3D-BTE calculations, however, the computational cost scales as $\mathcal{O}(n_{\rm unitcell}^{4}M^{2})$, where $n_{\rm unitcell}$ is the number of atoms in the unit cell and $M$ is the number of wave vectors in the first Brillouin zone.Turney et al. (2009) As a result, although we can perform MD simulations for films with thicknesses up to 100 unit cells for argon (25,600 atoms) and 40 unit cells for silicon (20,480 atoms), the 2D-BTE calculations are limited to thicknesses of 20 unit cells for argon and 8 unit cells for silicon.

II.2 Molecular Dynamics Simulation

The MD simulations on the bulk and film structures are performed using the open-source LAMMPSPlimpton (1995) and GPUMDFan et al. (2017) packages. LAMMPS is applied to obtain the finite temperature atomic structures and the RMS displacements for both argon and silicon and to predict the thermal conductivity of argon using the Green-Kubo method. For the silicon thermal conductivity predictions, GPUMD is used to calculate the heat current vector required in the Green-Kubo method because it provides the correct formulation for the Tersoff three-body potential.Fan et al. (2015) The time-step is 2 fs for both argon and silicon.

To obtain the zero-pressure bulk atomic structures, an MD simulation is run on an $8\times 8\times 8$ system in the $NVT$ ensemble (a constant number of atoms, volume, and temperature) for $5\times 10^{5}$ time steps to set the temperature. The system is then evolved in the $NPT$ ensemble (a constant number of atoms, pressure, and temperature) at zero pressure and the set temperature for $2\times 10^{6}$ time steps. The zero-pressure lattice constant is obtained by time-averaging the simulation cell size in the $x$-, $y$-, and $z$- directions and dividing by $N_{x}=N_{y}=N_{z}=8$. The simulation cell size is collected every 10 time steps and the averaging is performed over the last $10^{6}$ time steps in the $NPT$ ensemble.

For each film thickness and temperature, the zero-pressure in-plane lattice constant is determined using the same procedure as for bulk, but by time-averaging the simulation cell size in the $x$- and $y$- directions. Fixing the in-plane lattice constant at this value, an $NVT$ simulation is run for $2\times 10^{6}$ time steps to allow the film to relax in the cross-plane direction. Every ten time steps, the average $z$-location of each atomic layer is obtained by averaging the $z$-positions of all the atoms in that layer. The results are then time-averaged over the last $10^{6}$ time steps. The resulting atomic structure is used as the starting point for all further MD simulations and the 2D-BTE calculations. The film thickness is defined as the product of the number of atomic layers and the average layer separation in the cross-plane direction.

Thermal conductivity is predicted using the Green-Kubo method, which is an equilibrium MD technique based on the fluctuation-dissipation theorem.McQuarrie (2000) The system is first equilibrated in the $NVT$ ensemble for $2\times 10^{6}$ time steps and then in the $NVE$ ensemble for $2.5\times 10^{6}$ time steps. The heat current vector is then recorded every five time steps in the $NVE$ ensemble for $2.5\times 10^{7}$ time steps for argon and $2.5\times 10^{8}$ time steps for silicon and then used to obtain the heat current autocorrelation function (HCACF). The correlation time is $2.5\times 10^{5}$ time steps for argon and $2.5\times 10^{6}$ time steps for silicon. The HCACFs for thirty seeds with randomized initial velocities are collected to sufficiently explore the phase space. The thermal conductivity is obtained by integrating the HCACF for each seed and then averaging over the last 25% of the correlation time. We determined the uncertainty by calculating the 95% confidence interval, $1.96\frac{\sigma}{\sqrt{n}}$, where $\sigma$ and $n$ are the standard deviation and number of seeds. The maximum uncertainty for both material systems for all film thicknesses and temperatures is 7%. In terms of size effects, an in-plane size of $N_{x}=N_{y}=8$ provides a thermal conductivity that is converged to within 5%.

II.3 Lattice Dynamics Calculations

All lattice dynamics and BTE calculations are performed using in-house codes.McGaughey et al. (2019); Jain (2015) From the linearized BTE and the Fourier law, the thermal conductivity in direction $\alpha$, $k_{\alpha}$, can be expressed as

[TABLE]

where a phonon mode is identified by its wave vector $\boldsymbol{\kappa}$ and a polarization branch $\nu$, $c$ is the volumetric specific heat, and $v_{\rm g,\alpha}$ is the $\alpha$-component of the group velocity vector, which can be obtained from the phonon dispersion. The phonon dispersion (i.e., the relationship between the frequencies and wave vectors) is calculated from harmonic lattice dynamics calculations under the quasi-harmonic approximation (i.e., using the finite-temperature structures, as described in Sec. II.2). The phonon lifetime, $\tau_{\alpha}$, corresponds to a heat flow in the $\alpha$-direction. Scattering due to three-phonon interactions is calculated using anharmonic lattice dynamics calculations combined with an iterative solution to the linearized BTE.Wallace (1972); Reissland (1973); Ward (2009); Jain (2015); McGaughey et al. (2019) Classical (i.e. Maxwell-Boltzmann) statistics are used to allow for comparison with the inherently classical MD simulations, such that $c_{\rm v}=k_{\mathrm{B}}/V$, where $k_{\mathrm{B}}$ is the Boltzmann constants and $V$ is the simulation cell volume. We do not consider any additional scattering mechanisms for the bulk systems or for films modeled using the 2D Brillouin zone (i.e., a unit cell that spans the film thickness, as shown in Fig. 2).

For films modeled using the bulk 3D Brillouin zone, we account for phonon-boundary scattering using the Fuchs-Sondheimer model. The lifetime is calculated by modifying the intrinsic lifetime $\tau_{\rm pp,\alpha}$ withTurney et al. (2010)

[TABLE]

where $F$ is a mode-dependent scaling factor (i.e., a suppression function) given by

[TABLE]

Here, $p\mspace{-4.0mu}\left(\mspace{-8.0mu}\begin{smallmatrix}&\mathbf{\boldsymbol{\kappa}}\\ &\nu\end{smallmatrix}\mspace{-3.0mu}\right)$ is the mode-dependent specularity parameter, which ranges from zero for completely diffuse scattering to unity for completely specular scattering. Unless noted, we set $p\mspace{-4.0mu}\left(\mspace{-8.0mu}\begin{smallmatrix}&\mathbf{\boldsymbol{\kappa}}\\ &\nu\end{smallmatrix}\mspace{-3.0mu}\right)$ to be zero. $\delta\mspace{-4.0mu}\left(\mspace{-8.0mu}\begin{smallmatrix}&\mathbf{\boldsymbol{\kappa}}\\ &\nu\end{smallmatrix}\mspace{-3.0mu}\right)$ is equal to $L_{z}/(|v_{{\rm g},z}|\tau_{\rm pp,\alpha})$ where $v_{{\rm g},z}$ is the cross-plane ($z$) component of the bulk group velocity vector and $L_{z}$ is the film thickness. The scaling factor $F\mspace{-4.0mu}\left(\mspace{-8.0mu}\begin{smallmatrix}&\mathbf{\boldsymbol{\kappa}}\\ &\nu\end{smallmatrix}\mspace{-3.0mu}\right)$ decreases monotonically with decreasing film thickness, thus leading to a lower thermal conductivity.

The harmonic (second-order) and cubic (third-order) force constants required in the 2D/3D-BTE calculations are obtained by applying displacements of 0.005 Å to selected atoms and performing finite differencing of the resulting forces.McGaughey et al. (2019) For argon, the harmonic (cubic) force constants are extracted by applying a cutoff radius to include up to the seventh (first) nearest-neighbors [1.02 nm (0.58 nm)]. For silicon, the harmonic cutoff radius is set as the potential cutoff (0.3 nm) and the anharmonic cutoff radius includes up to the third nearest-neighbor (0.43 nm). Phonon properties and thermal conductivity are calculated on a $24\times 24\times 24$ wave vector grid for bulk and a $24\times 24$ wave vector grid for the films. Further increasing the resolution changes the predicted thermal conductivities by less than 1%.

III Lennard-Jones Argon Results

III.1 Film Atomic Structure

The in-plane lattice constants of the argon films with thicknesses between two and ten unit cells (1.1 - 5.3 nm) at temperatures of 5, 10, and 20 K are plotted in Fig. 3(a). Also plotted are the corresponding bulk lattice constants, which are provided in Table 1. At all temperatures, the in-plane film lattice constant increases monotonically with increasing thickness and approaches the bulk value.

The cross-plane atomic layer separations for films with thicknesses between two and ten unit cells at a temperature of 5 K are plotted in Fig. 3(b). Note that one unit cell is made up of two layers. Also plotted is the bulk layer separation, which corresponds to half of the bulk lattice constant. The layer separation varies across the film thickness and, as expected, is symmetrical about the centerline. The outermost layer separation for all films is larger than that in middle due to the weak bonding of the LJ potential. The layer separations inside the films are close to the bulk value. The exception is the two unit cell film, where the two middle layers are directly adjacent to the surface layers. The deviation from the bulk layer separation for this film at the centerline is 1%, which is large compared to the 0.04% deviation for the 10 unit cell film. Due to the cross-plane expansion and the in-plane contraction, the resulting film structures are anisotropic, an effect that becomes stronger as the thickness decreases. Similar results were found at temperatures of 10 and 20 K, as shown in Figs. S1(a) and S1(b) of the Supplemental Material.SM The anisotropy suggests that the anharmonicity present in the films will be different from that in bulk, which we explore next in Sec. III.2.

III.2 Anharmonicity

Phonon-phonon scattering is a result of anharmonicity and is the origin of the finite thermal conductivities of our bulk structures modeled using MD-GK and 3D-BTE and our film structures modeled using 2D-BTE. Before presenting the thermal conductivity predictions in Sec. III.3, we first examine the anharmonicity by considering the atomic RMS displacements.Parrish et al. (2014) As an atom deviates further from its equilibrium position, it is exposed to greater anharmonicity. Increasing anharmonicity is typically associated with increasing temperature, but it can also be a result of a free surface.

The RMS displacements are calculated from MD simulations in the $NVT$ ensemble. Data were collected every ten time steps from a $2\times 10^{6}$ time step simulation and averaged over each layer. The in-plane and cross-plane RMS displacements for film thicknesses between two and ten unit cells at a temperature of 5 K are plotted in Fig. 4. The in-plane RMS displacement is averaged over the $x$- and $y$- directions. Also plotted is the bulk value along one of the Cartesian directions, which is provided in Table 1. The in-plane and cross-plane RMS displacements increase in moving from the centerline to the surfaces. At the surfaces, the cross-plane RMS displacement is larger than the in-plane value, which is a result of the structural anisotropy discussed in Sec. III.1.

From Fig. 3(b), for film thicknesses greater than two unit cells, the two outermost layers have larger-than-bulk separations while the other layers are bulk-like. For the RMS displacements, we see that the effect of the surfaces penetrates three layers deep. The two unit cell film has the largest RMS displacements as it contains no bulk-like atoms. The larger displacements of the surface atoms will lead to larger anharmonicity in that part of the film. As the film thickness increases above five unit cells, the surface atom RMS displacements converge, indicated by the dashed and dotted lines for cross-plane and in-plane in Fig. 4. The converged in-plane and cross-plane values are 0.0065 nm and 0.0074 nm, which are 35% and 54% larger than the bulk value. Similar results are observed at temperatures of 10 and 20 K, as shown in Figs. S2(a) and S2(b).

III.3 Thermal Conductivity

The thermal conductivities of bulk argon predicted from the MD-GK and 3D-BTE methods at temperatures of 5, 10, and 20 K are provided in Table 1. The two methods can be compared based on how they treat the anharmonicity. MD simulations include the full anharmonicity and thus all orders of phonon-phonon interactions. We thus take the MD-GK thermal conductivities as a benchmark. The 3D-BTE and 2D-BTE approaches, as implemented, only consider up to three-phonon processes. When considering equivalent structures, potential parameters, and temperatures, the thermal conductivity prediction from 3D-BTE or 2D-BTE should therefore be larger than that from MD-GK, with a difference that increases with increasing temperature. The data in Table 1 follow this trend. Of note is the agreement at a temperature of 5 K, indicating that three-phonon processes are sufficient to capture phonon transport in the bulk phase.

The in-plane thermal conductivity variation with film thickness for argon films predicted using the 2D-BTE and MD-GK methods at a temperature of 5 K is plotted in Fig. 5(a). The bulk thermal conductivities from Table 1 are included as horizontal lines. Both methods predict a decrease in thermal conductivity with decreasing thickness up to a thickness of four unit cells ($2.1$ nm). Below that thickness, the thermal conductivity plateaus. The origin of the plateau is investigated in Sec. III.4 by analyzing the mode-level phonon properties. The 2D-BTE film thermal conductivities are always higher than the MD-GK values, with a difference that increases with decreasing thickness. The increase of the difference is a result of the increased anharmonicity at the film surfaces, which becomes more prominent as the thickness decreases, as discussed in Sec. III.2. The 2D-BTE calculations, which are based on small deviations from the equilibrium positions, do not capture the increased anharmonicity and thus over-predict the thermal conductivity for the smallest thicknesses. In using the MD-GK method to model LJ argon films at a temperature of 20 K, Turney et al.Turney et al. (2010) observed a monotonically decreasing thermal conductivity down to thicknesses as small as two unit cells. One possible origin of this different trend compared to our result is that Turney et al. used a fixed in-plane lattice constant equal to the bulk value,rep while we relaxed the film structures to zero-pressure as described in Sec. II.2.

Film thermal conductivities predicted from 3D-BTE combined with the Fuchs-Sondheimer phonon-boundary scattering model are also plotted in Fig. 5(a). Three data sets are shown. In the first data set, denoted as 3D-BTE-FS ($p=0$) and plotted as inverted brown triangles, the full set of phonon modes from the Brillouin zone are used and the phonon-boundary scattering is taken as completely diffuse. For the two unit cell film, the thermal conductivity from 3D-BTE-FS ($p=0$) is 68% smaller than that from the rigorous 2D-BTE treatment. The underestimation decreases as the film thickness increases and its modes become more bulk-like. In the second data set, denoted as 3D-BTE-FS ($p=0.62$) and plotted as inverted purple triangles, the specularity parameter $p$ is increased to 0.62 to match the 2D-BTE data for thicknesses greater than ten unit cells. With this choice, the plateau at very small thicknesses is not observed. We note that while the specularity parameter is in general phonon mode-dependent, we consider a mode-independent value to examine its overall impact on the predictions.

In the third data set, denoted as Depletion+FS and plotted as filled grey circles, we consider the phonon depletion approach discussed in Sec. ITurney et al. (2010); Wang and Huang (2014) with completely diffuse surfaces. The phonon wave vectors in the cross-plane direction are limited to discrete values by the constraint

[TABLE]

where $l_{z}$ is an integer with a magnitude less than $N_{z}/2$ and $a$ is the bulk lattice constant. The phonon properties are then calculated using the reduced list of wave vectors and the Fuchs-Sondheimer model. For the largest film considered using Depletion+FS ($N_{z}=24$), the predicted thermal conductivity matches the 3D-BTE-FS ($p=0$) result due to a Brillouin zone resolution that is the same as that used in the bulk calculation. As the film thickness is decreased, the thermal conductivity initially stays constant and then increases. It reaches a value close to the bulk thermal conductivity when $N_{z}=4$, which is clearly a non-physical result. The trend of increasing thermal conductivity with decreasing film thickness for the Depletion+FS calculation is consistent with the results of Turney et al., who studied LJ argon films at a temperature of 20 K.Turney et al. (2010) Although phonon depletion is an intuitive idea and a convenient way to consider phonon confinement, it remains an approximate 3D treatment that is not able to capture the correct trends or magnitudes of the thermal conductivity of ultrathin films.

The MD-GK, 2D-BTE, and 3D-BTE-FS thermal conductivity prediction methods were also applied at temperatures of 10 and 20 K. The results are plotted in Figs. S3 and 5(b). For MD-GK and 2D-BTE at both temperatures, the thermal conductivity plateau at very small thickness is again present, with the 2D-BTE results higher than those from MD-GK. Beyond a thickness of 10 unit cells at a temperature of 20 K, the MD-GK thermal conductivities unexpectedly become greater than the 2D-BTE values. This crossover, which we attribute to the use of the MD-derived structure in the 2D-BTE calculations, is discussed in Sec. S3. To match the 2D-BTE predictions at a temperature of 10 K, the 3D-BTE-FS calculations require a specularity parameter of 0.34. At a temperature of 20 K, completely diffuse boundaries are sufficient to provide a match between the 2D-BTE and 3D-BTE-FS predictions. As temperature increases, the atomic RMS displacements increase (Fig. 4), making the surface rougher and decreasing the specularity parameter. The decrease of the fitted specularity parameter with increasing temperature follows this trend. Furthermore, the impact of the specularity parameter decreases with increasing temperature, where shorter intrinsic MFPs limit the reduction of the effective MFPs due to phonon-boundary scattering.

III.4 Signatures of Confinement

We now analyze the phonon mode properties to identify signatures of confinement, with a focus on the 5 K system. The [100] phonon dispersion curves for the bulk and two unit cell film systems are plotted in Fig. S4. The phonon DOS for the bulk, two unit cell film, and fifteen unit cell film systems are plotted in Fig. 6. The phonon dispersion and DOS were obtained by applying harmonic lattice dynamics calculations to the the time-averaged MD structures described in Sec. III.1. The film dispersions are characterized by the emergence of a quadratic acoustic flexural mode, as found in 2D materials. The bulk and fifteen unit cell film DOS are indistinguishable, with a maximum frequency of 12.5 Trad/s. The similarities in these two DOS suggests that the thermal conductivity reduction for the fifteen unit cell film is due to changes in the phonon scattering. For the two unit cell film, however, the maximum frequency shifts downwards to 11.8 Trad/s and there is an increase in the DOS between frequencies of 3.0 and 5.2 Trad/s. Kim et al. link these changes to the film’s free surfaces, which remove the degeneracy of the transverse dispersion branches along high-symmetry directions.Kim et al. (2019) The peak in the DOS at a frequency of 4.9 Trad/s in the two unit cell film corresponds to a flattened phonon branch, as shown in Fig. S4.

We label phonons modes with frequencies below 6 Trad/s as low-frequency. This cutoff is chosen because it corresponds to the end of the quadratic portion of the bulk DOS. Low-frequency phonon modes in the two unit cell film account for 37% of the total population. As the thickness increases, the percentage of low-frequency phonon modes decreases and approaches the bulk value of 19%. We hypothesize that the large DOS at low frequencies in the very thin films is a signature of confinement.

To test this hypothesis, the thermal conductivity accumulations as a function of frequency for bulk (3D-BTE) and for films (2D-BTE) with thicknesses between two and fifteen unit cells are plotted in Fig. 7(a). The vertical coordinate of any point on the accumulation function represents the thermal conductivity contribution from phonon modes with frequencies less than the horizontal coordinate of that point (i.e., it is a cumulative distribution function). The vertical dashed line denotes the cutoff for the low-frequency phonon modes.

From bulk to films, the accumulation function at all frequencies decreases as the thickness decreases, particularly the contribution from high-frequency phonons, until the thickness reaches 10 unit cells. This trend is what would be expected based on a boundary scattering model like Fuchs-Sondheimer. Below this point, the accumulation functions in the low-frequency region increase with decreasing thickness. For the two-, three, and four unit cell films, the shape of the accumulation function is significantly different from that for bulk.

To further analyze this behavior, the contribution of the low-frequency phonon modes to thermal conductivity is plotted in Fig. 7(b) as a function of film thickness. For bulk, phonon modes with frequencies below 6 Trad/s (i.e., what we have identified as low-frequency) contribute 2.5 W/m-K to the total thermal conductivity of 9.0 W/m-K, which is plotted as a solid horizontal line. For the fifteen unit cell film, the contribution is 1.7 W/m-K, which is lower than bulk due to phonon-boundary scattering. Below a thickness of ten unit cells, however, the thermal conductivity contribution from the low-frequency phonons increases as the film thickness decreases. The contribution reaches a value of 2.6 W/m-K for the two unit cell film, which is larger than that for bulk. This behavior cannot be predicted from a model that includes phonon-boundary scattering and is a clear indication of confinement. The increased contribution of low-frequency phonon modes directly correlates with the increased DOS of these modes, as shown in Fig. 6. Interestingly, as also shown in Fig. 7(b), the contribution of low-frequency phonon modes to the thermal conductivities of argon films at temperatures of 10 and 20 K also starts to increase below a thickness of ten unit cells. Recall that all of our lattice dynamics-based predictions are made using classical statistics so that they can be compared to the predictions from the MD simulations. If quantum statistics are used, such that low-frequency modes are more populated than high-frequency modes, the effect of confinement is even stronger, as shown in Fig. S5.

We conclude that, for argon, independent of temperature: (i) When the film thickness is less than ten unit cells, phonon confinement leads to an increased DOS and accumulated thermal conductivity for low-frequency modes. When this increase is sufficiently large, the total thermal conductivity plateaus. (ii) When the thickness is greater than ten unit cells, the confinement effect is limited and phonon-boundary scattering is dominant.

IV Silicon Results

IV.1 Film Atomic Structure

The in-plane lattice constant for the silicon films is plotted in Fig. 8(a) for thicknesses between two and ten unit cells (1.0 - 5.5 nm). The bulk lattice constant (0.5440 nm) is plotted as a horizontal line. The films expand in the in-plane direction as the thickness is reduced, which is opposite to the argon result shown in Fig. 3(a). The cross-plane layer separations are plotted in Fig. 8(b). The bulk layer separation is plotted as a horizontal line. For thicker films, similar to the argon results shown in Fig. 3(b), the layers around the centerline are bulk-like. As the thickness decreases, the deviation from bulk at the centerline increases to 1.5% for the two unit cell film, larger than that of 1% observed in Fig. 3(b) for argon. At the surfaces, the layer separation is smaller than the bulk value, which is opposite to the findings for argon. Only the outermost layer shows a contraction, while the second outermost layer has a separation that is 0.2% larger than bulk. The contraction at the surface is a result of strong attractive forces from the first few neighbor shells, which originate from the stiffness of the Tersoff potential. The silicon and argon films are thus both anisotropic, but with changes in opposite directions compared to the bulk.

IV.2 Anharmonicity

The cross-plane atomic RMS displacements for silicon films with thicknesses between two and six unit cells are plotted versus layer number in Fig. 9 along with the bulk value. The in-plane values are provided in Fig. S6. Similar to argon, both surface atoms and interior atoms in films with small thickness have larger RMS displacements than in bulk, which is an indication of stronger anharmonicity. The RMS displacements of the surface atoms are larger than those of the interior atoms, even though their layer separation is smaller [Fig. 8(b)]. As the thickness increases, the RMS displacement of the atoms in the central interior layer approach the bulk value of 0.0054 nm, while that for the surface atoms reach 0.0088 nm (shown as a dashed horizontal line), which is 67% larger. The expansion of the surface layers in argon and their contraction in silicon both increase anharmonicity. This observation is consistent with the findings of Parrish et al., who found an increase in anharmonicity induced by tension in bulk argon and by compression in bulk silicon.Parrish et al. (2014)

IV.3 Thermal Conductivity

The thermal conductivities of the silicon films predicted from the MD-GK, 2D-BTE, and 3D-BTE-FS methods are plotted in Fig. 10. To allow for a better comparison with Fig. 1, the thermal conductivities are scaled by their corresponding bulk values: 414 W/m-K for MD-GK and 464 W/m-K for 3D-BTE. The bulk predictions are higher than the experimental value of 293 W/m-K.Inyushkin et al. (2004) The Tersoff potential is well-known for overpredicting thermal conductivity.Broido et al. (2005) The raw thermal conductivities are plotted in Fig. S7.

As the thickness decreases, both MD-GK and 2D-BTE predict a decreasing thermal conductivity. There is a minimum at a thickness of three unit cells, below which the thermal conductivity shows an unexpected increase of about 15%. This behavior is similar to the argon data plotted in Fig. 5, but with an increase instead of a plateau. This result suggests that the phonon confinement effect is stronger in silicon, a hypothesis that we further explore in Sec. IV.4 by considering the mode-dependent phonon properties.

As expected, the 2D-BTE thermal conductivities are always higher than the MD-GK values. The deviation increases as the film thickness decreases and reaches a maximum of 91 W/m-K for the 2 unit cell film, as shown in Fig. S7. The increase in the deviation is analogous to that in argon as discussed in Sec. III.3, and is consistent with the increase of the anharmonicity as the film gets thinner, as shown in Fig. 9.

The scaled thermal conductivities predicted from 3D-BTE-FS ($p=0$) (i.e., with a specularity parameter of zero) vary from 0.08 (two unit cell film) to 0.22 (forty unit cell film), which is much lower than 2D-BTE scaled values, which range from 0.53 to 0.62 over the computationally accessible range of two to eight unit cells. Based on the argon results, the potential role of specular phonon-boundary interactions must be assessed. We first consider the experimentally-derived wavelength-dependent specularity parameters for silicon reported by Ravichandran et al. at a temperature of 300 K.Ravichandran et al. (2018) We use the upper bound of their results, which vary from 0.18 to 1 as the wavelength increases. Applying their results to the 3D-BTE-FS model increases thermal conductivity, leading to scaled values ranging from 0.12 (two unit cell film) to 0.30 (forty unit cell film). These values remain smaller than the exact 2D-BTE treatment. To achieve a good agreement between 3D-BTE-FS and 2D-BTE, a mode-independent specularity parameter of 0.98 is required, which is non-physical at finite temperatures for the full phonon spectrum. The thermal conductivity increase at the smallest thickness is not observed with any choice of specularity parameter, consistent with the results for argon.

The scaled thermal conductivities predicted from MD-GK range from 0.36 [the minimum value at three unit cells (1.6 nm)] to 0.76 [for the maximum thickness considered of forty unit cells (21.8 nm)]. Neogi et al.Neogi et al. (2015) applied MD-GK to Tersoff silicon films with reconstructed $2\times 1$ surfaces. They predicted scaled thermal conductivities at a temperature of 300 K that varied from 0.37 (1 nm) to 0.71 (20 nm), consistent with our results. Their MD-GK predictions for these silicon films were ten times larger than their experimental measurements (green triangles in Fig. 1) and those of Chávez-Ángel et al. (black triangles in Fig. 1). Chávez-Ángel et al. (2014) By including surface roughness and a native oxide layer in their MD simulations, they obtained scaled thermal conductivities that agreed with the experimental measurements. They then etched the experimental oxide layer and measured increased thermal conductivities (blue triangles in Fig. 1), which agreed with their MD simulations of rough films without the oxide. As reported by Xiong et al.,Xiong et al. (2017) the native oxide layer can introduce local resonant modes that interact with the propagating phonons inside the film, reducing their group velocity and lowering thermal conductivity. This phenomenon is analogous to predictions for nanophononic metamaterials.Honarvar and Hussein (2018); Davis and Hussein (2014)

As shown in Fig. 1, thermal conductivity predictions from 3D-BTE-FS ($p=0$) (corresponding to the FS+Esfarjani curve) agree with the experimental TTG data for suspended silicon films with thicknesses as small as 20 nm. 3D-BTE-FS does not include the effects of confinement, however, and, as we showed in Fig. 10, leads to a large underestimation of thermal conductivity compared to the rigorous 2D-BTE and MD-GK approaches, even with realistic specularity parameters. As found by Neogi et al.,Neogi et al. (2015) the low thermal conductivities of suspended silicon films measured by Raman thermometry with thicknesses smaller than $\sim 200$ nm (green and black triangles in Fig. 1) result from the effects of confinement and the presence of a native oxide layer. None of the Raman thermometry measurements are well predicted by 3D-BTE-FS, which does not capture thermal transport inside the oxide and/or the impact of the oxide on the film’s phonon properties. The same holds true for the suspended island measurements (red and pink squares). It is thus possible that the reasonable agreement found between 3D-BTE-FS ($p=0$) and the TTG experimental measurements for ultrathin silicon films is coincidental due to use of bulk phonon properties with a boundary scattering modeling and ignoring the oxide layer.

IV.4 Signatures of Confinement

Following the argon analysis presented in Sec. III.4, the DOS of bulk and silicon film systems (two and six unit cells) are plotted in Fig. 11. The quadratic portion of bulk DOS ends at a frequency of 31 Trad/s, which we identify as the low-frequency region. There are two peaks between frequencies of 20 and 25 Trad/s in the films that are not present in bulk, similar to argon. For the two unit cell film, 23% of the phonon modes are low-frequency, which is more than double that for bulk.

The frequency-dependent thermal conductivity accumulation functions for silicon bulk and films are plotted in Fig. 12(a). In both bulk and films, phonon modes with low frequencies make the dominant contribution to thermal transport. In moving from bulk to the films, the accumulation function decreases at all frequencies as the thickness decreases, until the thickness reaches six unit cells. Below this point, the contributions from phonon modes with frequencies below 31 Trad/s (i.e. low-frequency phonon modes) increase as the thickness is further decreased. This trend is similar to that observed in argon. Furthermore, the accumulation function for the two unit cell film is significantly different compared to bulk.

The contribution of low-frequency phonon modes to thermal conductivity is plotted in Fig. 12(b). For bulk, the low-frequency phonon modes contribute 222 W/m-K to the total thermal conductivity of 464 W/m-K. For the eight unit cell film, the contribution is 157 W/m-K, which is lower than bulk due to phonon-boundary scattering. Below a thickness of six unit cells, the contribution from the low-frequency phonon modes increases as the film thickness decreases and they exhibit a larger DOS. The contribution reaches a value of 242 W/m-K for the two unit cell film, which is larger than that for bulk. This trend is similar to that observed in argon [Fig.7(b)]. Interestingly, the transition point in silicon is 24 atomic layers, which is close to the 20 atomic layers (10 unit cells) observed for argon.

Low-frequency phonon modes account for 48% of the bulk silicon thermal conductivity and 28% of the bulk argon thermal conductivity. For the two unit cell films, the contribution of the low-frequency phonon modes is 85% for silicon and 67% for argon. We believe that the greater contribution for the silicon films is what leads to the increase in their total thermal conductivity for the smallest thicknesses, compared to the plateau observed for argon.

V Summary

We investigated in-plane thermal transport by phonons in ultrathin films of (i) LJ argon at temperatures of 5, 10, and 20 K, and (ii) Tersoff silicon at a temperature of 200 K using MD simulations, lattice dynamics calculations based on 2D and 3D Brillouin zones, and the BTE. The finite-temperature atomic structures used for all thermal conductivity predictions were generated from MD simulations. As shown in Figs. 3 and 8, there is anisotropy in the in-plane and cross-plane atomic spacings in the films compared to the corresponding isotropic bulk structures. RMS atomic displacement obtained from MD simulations revealed that surface atoms experience a stronger anharmonicity than those at the film center, as shown in Figs. 4 and 9, which has important consequences when comparing the MD-GK and 2D-BTE thermal conductivity predictions.

For argon modeled with MD and 2D lattice dynamics, a thermal conductivity plateau is observed for films with thicknesses smaller than four unit cells (2.1 nm) (Fig. 5). For the same calculations for silicon, the film thermal conductivity increases when its thickness is reduced below three unit cells (1.6 nm) (Fig. 10). We determined that these unexpected trends are a signature of phonon confinement, which causes an increase in the DOS of low-frequency phonons and a corresponding increase in their contribution to thermal conductivity to a level that can exceed that present in the bulk structures, as shown in Figs. 7(b) and 12(b). This effect could potentially be observed experimentally by applying emerging techniques for measuring phonon mode-dependent contributions to thermal conductivity.Regner et al. (2015)

Bulk phonon properties combined with the Fuchs-Sondheimer model for phonon-boundary scattering were also used to predict the film thermal conductivities (i.e., 3D-BTE-FS). These predictions did not show the plateau (for argon) or the increase (for silicon) in thermal conductivity at very small thicknesses. The common assumption of fully diffuse phonon-boundary scattering was not able to capture the thermal conductivity magnitudes for thicker argon films at temperatures of 5 and 10 K or for the silicon film at a temperature of 200 K. We also showed that the phonon depletion method produces non-physical thermal conductivities.

Our work demonstrates that the rigorous 2D description of phonons in ultrathin films is essential for accurate thermal conductivity prediction and for elucidating the underlying transport physics. As advances continue to be made for the fabrication and characterization of ultrathin films, and notably suspended films,Cuffe et al. (2015); Neogi et al. (2015); Chávez-Ángel et al. (2014); Wang et al. (2017) the tools described herein will thus be of utmost importance.

Acknowledgements.

B. F. and G. T. acknowledge financial support from the National Natural Science Foundation of China under grant number 51825604. K. D. P., H.-Y. K., and A. J. H. M. acknowledge financial support from the National Science Foundation Award DMR-1507325. B. F. acknowledges financial support from a program of the China Scholarships Council (Award No. 201706280251). G. T. acknowledges financial support from the 111 Project under grant number B16038.