{'fraction_non_alphanumeric': 0.0516388242757454, 'fraction_numerical': 0.0222034256713893, 'mean_word_length': 4.655584788328151, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Size Dependence and Ballistic Limits of Thermal Transport in Anisotropic Layered Two-Dimensional Materials

Layered materials have uncommonly anisotropic thermal properties due to their strong in-plane covalent bonds and weak out-of-plane van der Waals interactions. Here we examine heat flow in graphene (graphite), h-BN, MoS2, and WS2 monolayers and bulk films, from diffusive to ballistic limits. We determine the ballistic thermal conductance limit (Gball) both in-plane and out-of-plane, based on full phonon dispersions from first-principles calculations. An overall phonon mean free path ({\lambda}) is expressed in terms of Gball and the diffusive thermal conductivity, consistent with kinetic theory if proper averaging of phonon group velocity is used. We obtain a size-dependent thermal conductivity k(L) in agreement with available experiments, and find that k(L) only converges to>90% of the diffusive thermal conductivity for sample sizes L>16{\lambda}, which ranges from ~140 nm for MoS2 cross-plane to ~10 um for suspended graphene in-plane. These results provide a deeper understanding of microscopic thermal transport, revealing that device scales below which thermal size effects should be taken into account are generally larger than previously thought.

Stimulated by extensive studies of graphene [1][2][3], much interest now exists in the properties of other two-dimensional (2D) layered materials [4][5][6][7], such as hexagonal boron nitride (h-BN) and transition metal dichalcogenides (TMDs, e.g., MoS2, and WS2). Understanding their thermal properties is important for improving thermally-limited performance in integrated electronics and energy conversion devices [8][9][10][11]. The thermal properties of 2D layered materials are unusual and highly anisotropic, due to strong in-plane chemical bonding and weak interlayer van der Waals interactions [1,12]. For instance, bulk graphite has over 300-fold higher in-plane than out-of-plane thermal conductivity at room temperature; in h-BN this anisotropy is 200-fold. Like in other materials, when sample dimensions become comparable to the average phonon scattering length or mean free path (MFP) λ, thermal transport can no longer be described with classical parameters (such as the bulk thermal conductivity) and so-called size effects become important [8].
Simply stated, if the sample size (i.e. "length") in the direction of heat flow L ≪ λ, heat conduction enters a ballistic regime (with no scattering) and the thermal conductance approaches a ballistic upper limit, Gball, which is only a function of temperature. Interestingly, if the classical relationship between conductance and conductivity is applied (k = GL/A) the thermal conductivity becomes sizedependent and approaches zero in the limit of the shortest ballistic samples. The ballistic thermal conductance of carbon nanotubes [13], silicon [14] and graphene in-plane [15] have been recently examined, but these limits are unknown in other 2D (layered) materials, and in their out-of-plane direction. More generally speaking, much less is known about the ballistic limits of thermal transport compared to those of electrical transport [16][17][18][19].
If the sample size L ≫ λ, heat conduction is said to be in the classical (Fourier) regime, with the thermal conductivity being independent of sample size. However, the transition from ballistic to diffusive is not well understood, nor is the apparent evolution of thermal conductivity with size, k(L), in this regime. Some theoretical models have predicted that in low-dimensional momentum-conserving systems k will diverge with L, scaling as ~L α with α = 0.22 to 0.5 for one-dimension (1D) [20][21][22] and as ~log(L) for 2D [22][23][24]. Although it is natural to assume that the transition from ballistic to diffusive should occur over sample sizes of the order of the phonon MFP, this value is often not estimated consistently in the literature. For example, the Si phonon MFP is ~40 nm at room temperature [25] based on a Debye heat capacity model with an average sound velocity, but ~300 nm when the phonon dispersion is taken into account [26]. Analysis of cumulative thermal conductivity as a function of MFP has suggested that phonons with MFP > 1 µm contribute ~40% to k in Si [27][28][29][30][31], so a "median" MFP of 0.5-1 µm has also been suggested [25,31]. Thus, an improved understanding of size effects and phonon MFP are crucial both in bulk and in atomically thin materials. While a range of MFPs provides insight into the spectral distribution of phonon contributions [27][28][29][30][31], providing a single MFP gives information in a concise way, being useful to evaluate size effects and transport in practical devices, where the device dimensions are known.
In this work, we study the ballistic-diffusive thermal conduction in monolayer graphene, h-BN, MoS2 and WS2, as well as their three-dimensional (3D) bulk counterparts. Based on full phonon dispersions obtained from ab initio simulations, we calculate the in-plane (for monolayer and bulk) and cross-plane (for bulk) ballistic thermal conductance Gball of these materials for the first time. Due to their stronger chemical bonding, graphene (graphite) and h-BN generally show higher Gball than MoS2 and WS2 above ~100 K. We then obtain L-dependent thermal conductivity of these materials by using a ballistic-diffusive model, which for suspended graphene shows good agreement with recent measurements [32] of k vs. L up to 9 µm. We also discuss how an estimate of the phonon MFP as a rigorous average of all phonon modes can be condensed to an expression including Gball and the diffusive thermal conductivity kdiff. We find that kinetic theory can reach the same MFP estimate as long as a correctly averaged phonon group velocity is used. Importantly for anisotropic layered materials, 2D and 1D forms of the kinetic theory should be used for in-plane and out-of-plane transport, respectively. Based on the ballistic-diffusive model, k converges to kdiff only around L ~ 100λ, at much larger sizes than commonly assumed. Thus, whether k is divergent in low-dimensional materials should be examined beyond ~100λ to distinguish from the intrinsic increase of k due to the ballistic-to-diffusive transition.

Phonon Dispersion.
To study the thermal properties of these systems, we begin by calculating the phonon dispersions using a first-principles approach. First, the equilibrium atomic structures are calculated within density functional theory (DFT) using projector-augmented wave (PAW) potentials [33], as implemented in the VASP code [34]. The local density approximation (LDA) is used for the exchange-correlation functional [35]. Then, the interatomic force constants are calculated in a supercell using the frozen-phonon approach [36]. Finally, phonon dispersions in the entire first Brillouin zone (BZ) are calculated by the PHONOPY package [37] interfaced with VASP. Details of these calculations are described in Supplemental Material Section 1.
The atomic lattices of bulk graphite, h-BN, and TMDs are shown in Figs. 1(a-c). Graphite tends to crystalize in the AB stacking, while h-BN favors the AA' stacking [38,39]. For TMDs, the most common 2H phase [5,40] (similar to the AA' stacking in h-BN) is considered here.
where ħ is the reduced Planck constant; q the sum with respect to q  is converted to an integral, and the ballistic thermal conductance per unit cross-sectional area (i.e., Gball/A with A  WH, where W and H are the material width and height, respectively) is calculated to give a value independent of size and easy for comparison among different materials. For 2D (monolayer) and 3D (bulk) layered materials, Gball/A is given, respectively, by where the integrals are over the entire 2D and 3D BZ, respectively. For monolayers, we take H = c/2, where c is the cross-plane lattice constant of bulk ( Fig. 1). For bulk, , ( ) v q  to give the in-plane (||) and cross-plane () Gball/A, respectively (see Supplemental Material Section 3 for the derivation and calculation details).
The calculated ballistic thermal conductance of monolayers, bulk in-plane (||) and cross-plane () as a function of temperature are shown in Fig. 4, for graphite, h-BN, MoS2 and WS2. For each material, the monolayer Gball/A is higher than that of bulk (||) at low T, but they become almost overlapped at high T, because the in-plane phonon dispersions of monolayer and bulk are nearly identical except at low frequencies. For bulk (), its Gball/A is similar to that of bulk (||) at very low T (< 10 K), but the in-plane (||) ballistic thermal conductance becomes ~10 times greater at all higher temperatures due to the significant anisotropy of the phonon dispersion.
To put these results in perspective, we also plot the ballistic thermal conductance of silicon, the most widely used semiconductor [14], as a dashed line in each panel of Fig The room-temperature Gball/A values and Debye temperatures ΘD [45][46][47][48] of these materials are listed in Table 1.
In the high T limit, all Gball/A vs. T curves flatten out (saturate) when all phonon modes are fully excited. The temperatures at which Gball/A reach 90% of their maximum value are nearly the same for the monolayer and bulk (||) of each material. These are ~1150 K, ~1060 K, ~320 K and ~310 K for graphite, h-BN, MoS2 and WS2, respectively, which are roughly 60% of their Debye temperature ΘD (Table 1). The corresponding temperatures for cross-plane bulk () are ~780 K, ~1270 K, ~345 K and ~360 K for graphite, h-BN, MoS2 and WS2, respectively. The same observation holds approximately for bulk Si as well, where 90% of its ballistic conductance is reached at ~400 K, approximately 60% of its ΘD (Table 1).
Interestingly, there is a "bump" in the cross-plane vz, unlike the nearly zero vz in other bulk materials. Thus, as temperature increases and these phonon modes start to contribute to cross-plane conduction, a sharp increase (bump) appears in the Gball/A vs.
T curve.
At low temperature, all Gball/A vs. T curves show power law scaling, ~T n (Fig. 4). For the monolayers, the power exponent is n ≈ 1.6−1.7, which is a combined effect of n = 1.5 from the quadratic ZA branch and n = 2 from the linear TA and LA branches [1,13,15]. The power law only applies to monolayers below ~90 K for graphene and h-BN, and below ~40 K for MoS2 and WS2. For bulk (||), the power exponent is n ≈ 2.9, which is a combined effect of n = 2.5 from the ZA branch (if purely quadratic) [13] and n = 3 from linear TA and LA branches (see Supplemental Material Section 8).
The reason that the overall n is closer to 3 than 2.5 is that the real ZA branch for bulk has a linear component (not purely quadratic) at very low frequency [49], ultimately leading to n > 2.5. The power law for bulk in-plane (||) is valid only below ~35 K for graphite and h-BN, and below ~20 K for MoS2 and WS2. For bulk cross-plane (), the power exponent is n ≈ 2.7−2.9 with the same physical reason as bulk in-plane (||), and it applies only below ~10 K for graphite and h-BN, and below ~7 K for MoS2 and WS2.
The calculated ballistic thermal conductance is a useful quantity in the investigation of heat conduction limits in nanomaterials and nanoscale devices. It gives the upper limit of heat flow at a given temperature, and also reflects the intrinsic anisotropy of heat conduction in such layered materials.
Any measured or calculated thermal conductivity (k) can be connected to conductance (G) as G/A = k/L, and then compared to Gball/A. The quantity k/L must always be less than or equal to Gball/A; if this relationship is not satisfied (as for some classical molecular dynamics simulations) the thermal conductivity is unphysical [50]. To put it another way, kball = GballL/A represents the maximum thermal conductivity a nanoscale device or sample can reach. In addition, the ratio of (k/L) to (Gball/A) gives the percentage of ballistic conduction that a nanoscale device reaches [15,32]. This approach can also be used to predict thermal conductivity as a function of device length L, and to estimate an overall phonon MFP λ averaged across all modes.
Length-Dependent Thermal Conductivity. We next turn to a discussion of length-dependent thermal conductivity in nanoscale layered materials. In the ballistic regime (L ≪ λ) the thermal conductance approaches a constant, thus the ballistic thermal conductivity kball = GballL/A becomes linearly dependent on length L. In the diffusive regime (L ≫ λ), the thermal conductivity is typically independent of length L (the case of k divergent with L will be discussed later), becoming a constant k = kdiff (diffusive thermal conductivity). In the diffusive regime the conductance scales inversely with sample length, as a thermal analogy to Ohm's law, G = kA/L.
Consequently, in the intermediate ballistic-diffusive (or quasi-ballistic) regime, the effective thermal conductivity should increase with L and gradually converge to kdiff; this behavior is similar to the effective mobility dependence on sample length during quasi-ballistic charge transport in shortchannel transistors [16,17]. This transition can be captured through a phenomenological ballistic-diffusive (BD) model [15,49]: where the first expression is a "1-color" model, and the second one is a branch-resolved "multicolor" model, taking into account different phonon branch (b) contributions. As we will show later, this is essentially a Landauer-like model and it can yield an expression of diffusive thermal conductivity consistent with kinetic theory.
We first apply this model to the thermal conductivity of graphene supported on SiO2. Having shown that the BD model can successfully describe the k change with L in supported graphene, next we focus on the L-dependent k in suspended graphene. Some theoretical studies have predicted that k will diverge with L in isolated low-dimensional (1D and 2D) momentum-conserving systems, i.e. k ~ L α for 1D [20][21][22] and k ~ log(L) for 2D like graphene [22][23][24]. However, these theoretical studies only consider pure 1D (2D) systems in which atoms can only move in one (two) dimensions, but atoms in realistic 1D (2D) materials like carbon nanotube (graphene) can still move in 3D, leading to properties that pure 1D and 2D systems cannot capture, such as radial breathing mode and flexural mode. It is hence unclear if predicted divergence will happen in realistic materials.
In addition, other studies [54][55][56] have argued that disorder and higher-order three-phonon scattering may eliminate the divergence and, as we will show, no experiments have conclusively observed divergent k yet, in neither 1D nor 2D systems. For supported graphene discussed above, the ~log(L) divergence does not appear due to phonon scattering with substrate vibrational modes [57,58], but for Assuming the divergence of k ~ log(L) for large L, the BD model can be adjusted with a logarithmic term to capture this effect (see Supplemental Material Section 8 for more discussion): where k0 and L0 are parameters with units of thermal conductivity and length, respectively. This modified "log" model reproduces the ballistic limit correctly as well, i.e. k → (Gball/A)L when L → 0.
Choosing kdiff = 1590 Wm -1 K -1 , k0 = 130 Wm -1 K -1 , and L0 = 1.1 µm, the log model can also fit the data of Xu et al. [32] within the bounds of experimental error [dash-dotted line in Fig. 5(b)]. Both models obey the ballistic limit and fit the available data, but they are distinguishable only at L > 10 µm, whereas data are available only up to ~9 µm. In the regime below 10 µm the increase of k with L mainly results from the ballistic-to-diffusive transition. This suggests that the effects of divergent thermal conductivity in suspended graphene will only show up unambiguously in samples longer than ~10 µm, and that suspended devices and samples shorter than this length can be adequately described with a simple, 1-color BD model.
A more general expression of k including both length and temperature dependence for SWCNTs is provided in our previous study of short channel (L < 100 nm) SWCNT devices [63], where heat conduction is nearly ballistic. MoS2 (1L) [68], MoS2 (||) [69], MoS2 () [70], and Si [71], respectively (see Table 1). These kdiff values are at room temperature, for undoped and isotopically unmodified samples. The layered materials show strong anisotropy in thermal conductivity, and their values span a wide range, more than three orders of magnitude. The estimated k dependence on L given here will be helpful for understanding heat conduction at the length scales where diffusive-ballistic effects take place. We note that convergence to "bulk" thermal conductivity does not occur until 10s of nanometers in the  direction and 100s of nanometers to microns in the || direction (up to 10 µm in suspended graphene).

Phonon Mean Free Path and Group
The transmission coefficient is λbs/(L + λbs), given in terms of the transport length L and an averaged phonon back-scattering MFP, λbs. A simple rearrangement of the first expression will reproduce the BD model, meanwhile yielding a relation λbs = kdiff/(Gball/A). This relation can be directly derived from the definition of a phonon MFP with all modes weighted properly across the phonon dispersion (see Supplemental Material Section 5). The common phonon MFP λ is shorter than the back-scattering MFP λbs, and they are related with a factor βλ (see Table 2), The factor βλ = 1/2, 2/π, and 3/4 for 1D, 2D, and 3D materials with isotropic thermal properties, respectively [73]. For anisotropic layered materials, the in-plane factor is βλ(||) = 2/π, same as isotropic 2D, while the cross-plane factor is βλ( ) = 1/2, same as 1D. The derivation is provided in Supplemental Material Section 5.
As long as reliable kdiff from calculations or experiments are available, the overall phonon MFP can be estimated based on Eq. (6), because Gball/A can be calculated without approximations. Taking graphite as an example, using its well-known measured kdiff(T) [1,65] and our calculated Gball(T)/A, the obtained phonon in-plane and cross-plane MFPs (λ|| and λ  ) as a function of temperature are shown as solid lines in Fig. 6(a). They decrease rapidly as T increases and differ by two and one orders of magnitude at low and high T, respectively. At room temperature, λ|| ≈ 290 nm and λ  ≈ 10 nm.
The latter corresponds to ~30-layers thickness in graphite, which means the cross-plane heat conduction in most few-layer graphene samples [74] is in the quasi-ballistic regime. The advantage of a single MFP averaged over all phonon modes [Eq. (6)] is giving information in a concise way, unlike the use of frequency-dependent MFPs [27][28][29][30][31], and it helps us understand when to make size-effect corrections to thermal conductivity in practical devices.
We also compare our estimates with the classical approach to evaluating the phonon MFP by kinetic theory as kdiff = (1/d)Cvvλ, where d = 1−3 is the material dimension, Cv is the heat capacity, and v is the phonon group velocity (typically taken as average sound velocity vs) [75]. The dashed lines in To obtain an improved phonon MFP by the classical kinetic theory, the group velocity v should be an average weighted by the heat capacity. By analyzing the integral expressions of Gball/A and Cv, we find the averaged group velocity can be simply given by where the velocity factor βv = 2, π, and 4 for isotropic 1D, 2D, and 3D materials. Similar to the MFP factor βλ, for anisotropic layered 3D materials the in-plane βv = π, same as for isotropic 2D, and the cross-plane βv = 2, same as for 1D [see Supplemental Material Section 7 for the derivation of Eq. (7) and βv]. Substituting Eq. (7) into the kinetic theory, the phonon MFP is then given by We can find d/βv ≡ βλ for all cases (see Table 2), thus when the correct v is used, the kinetic theory can give the same estimation of MFP as our model [Eq. (6)]. After demonstrating the BD model and kinetic theory can be reconciled in the MFP estimation, we apply Eq. (6) to other interesting materials whose measured kdiff are available.  Table 1.
By comparing the obtained λ with k(L) in Fig. 5(d), we can find that k reaches ~85% of kdiff at L ≈ 10λ and ~90% at 16λ. At larger device dimensions k increases slowly and reaches ~98% of kdiff (essentially fully diffusive) at L = 100λ, which is much longer than the commonly assumed sample sizes. The entire ballistic-to-diffusive transition takes more than two orders of magnitude in sample size (length), from L ≈ λ to 100λ to complete. This also indicates that any obtained k increase in this size regime is likely to arise from the ballistic-to-diffusive transition, and whether k is divergent in low-dimensional materials should be examined in samples with L > 100λ.
We note that Eq. (6) gives an MFP including contributions from all phonon modes, but some modes (e.g., optical) have smaller contributions to thermal conductivity and short MFPs (<10 nm for in-plane and <1 nm for cross-plane). This implies that the modes contributing significantly to heat conduction should have MFPs longer than the value given by Eq. (6). Indeed, by analyzing the cumulative thermal conductivity as a function of MFP, it was found that phonons with MFPs longer than 1 µm and 100 nm contribute ~40% to the total thermal conductivity for Si [27][28][29][30][31] and cross-plane graphite [76] at 300 K, respectively. These length scales correspond to a "median" MFP of ~10λ based on our definition from Eq. (6). Our approach provides a single, rigorous MFP with all phonon modes weighted properly (see Supplemental Material Section 5), which should be used in the kinetic theory and to understand quasi-ballistic thermal transport. From Fig. 5(d), when L < λ, heat conduction enters the ballistic regime (k linear with L) as expected. When L is shorter than the "median" MFP suggested above (~10λ), k is still a strong function of L, but is not in the ballistic regime, implying that this is a characteristic length scale below which the size effects (not just ballistic transport) take place.
In conclusion, we calculated the ballistic thermal conductance Gball of graphene and graphite, h-BN, MoS2, and WS2 based on full phonon dispersions obtained from ab initio simulations. We obtained the size-dependent thermal conductivity k with a simple ballistic-diffusive model verified against Boltzmann transport calculations. We showed that a single phonon MFP λ can be simply determined by Gball and kdiff, and that kinetic theory can reach the same result if a proper phonon group velocity is used. The MFP λ provides nanoscale heat flow information in a concise way, which is of great utility for simple and rapid estimates in practical nanoscale devices. Based on the calculated λ and our approach, we find that k converges to 90% kdiff at L ~ 16λ, and full convergence is only achieved for sample dimensions L > 100λ, much longer than commonly assumed. This also implies that to verify divergent k in low-dimensional systems, simulations and measurements should be performed for L > 100λ.   [41], red circles [42], and blue triangles [43] for graphite, and other experimental data shown for h-BN [44] and MoS2 [45].) Lines are the calculated phonon dispersions and symbols are available experimental data for bulk inplane modes. (Data for graphene are represented as red circles [42] and blue triangles [43]; experimental data are also shown for h-BN [44] and MoS2 [45].) The 2D Brillouin zone for monolayers is displayed as the inset in (a).