Reexamination of hydrodynamic phonon transport in thin graphite

The strong covalent bonding of light elements in graphite and graphene leads to high basal-plane lattice thermal conductivity.¹ Several recent theories have suggested that the high thermal conductivity consists of a large contribution from the low-frequency flexural phonons because of their large density of states and strong normal scattering.^2,3 In addition, the reflection symmetry in graphene prohibits phonon–phonon scattering processes that involve an odd number of flexural phonons. Normal scattering of the small-wavevector flexural phonons conserve the phonon momentum similar to intermolecular scattering. In an intermediate temperature range that is large compared to those found for other bulk crystals,^4–8 normal scattering of a large population of flexural phonons in the highly anisotropic layered structure is more frequent than Umklapp processes that do not conserve phonon momentum.^9,10 Consequently, a hydrodynamic phonon transport regime emerges in between the low-temperature ballistic regime without phonon–phonon scattering and the high-temperature diffusive regime dominated by momentum-destroying Umklapp or defect scattering.^9–12 The hydrodynamics theory of phonon transport in graphite has recently been supported by the observation of the second sound, which is the phonon equivalent of the sound wave in a gas, by a diffraction grating measurement¹³ and a pump-probe thermal reflectance measurement¹⁴ at temperatures as high as about 100 K.

In addition to second sound, hydrodynamic phonon transport can produce a similar heat flux pattern as the Poiseuille flow when the transport cross section is reduced to be comparable to the normal scattering mean free paths (MFPs).^15–17 Such phonon Poiseuille flow can give rise to a more rapid increase of the thermal conductance with temperature than the ballistic thermal conductance because momentum-conserving normal processes increase the momentum relaxation length compared to the ballistic limit. Similarly, it can lead to size-dependent thermal conductivity behaviors that are distinct compared to the size effects associated with ballistic transport. The size-dependent thermal conductivity caused by hydrodynamic, ballistic, and other non-diffusive phonon transport phenomena has received considerable interest in recent years due to both the rich fundamental physics and their relevance to heat dissipation in functional devices and materials.

Experimental studies of thermal transport in thin graphite and its monolayer graphene limit have continued to produce intriguing results. Electro-thermal measurements have observed decreased basal-plane thermal conductivity of few-layer graphene (FLG) in contact with a disordered support that scatters graphene phonons,^18,19 somewhat similar to conventional diffuse surface scattering that reduces the in-plane thermal conductivity of thin films and nanowires.²⁰ In comparison, a micro-Raman thermometry measurement obtained an increased basal-plane thermal conductivity as the layer thickness of suspended few-layer graphene (FLG) was decreased toward a single layer,²¹ in agreement with the expectation that elimination of phonon modes with out-of-plane wavevectors reduces the phonon scattering phase space in two-dimensional (2D) systems with defect-free, atomically smooth and specular surfaces.³ Notably, a recent thermocouple-based steady-state measurement reported an increase of both the thermal conductivity and the peak temperature of graphite with decreasing thickness from 580 to 8.5 μm.²² The observed trend is opposite to those found in other thin films with diffuse surfaces, which reduce the thermal conductivity magnitude and increase the peak temperature with decreasing thickness due to increased phonon–boundary scattering compared to the intrinsic phonon–phonon scattering processes.²⁰ Meanwhile, the thickness range for the reported thermal conductivity increase is orders of magnitude larger than that found in a prior first-principles phonon transport calculation,³ for which the calculated thermal conductivity of graphite starts to increase only when the thickness is reduced to about five atomic layers. Since the size confinement effects on the phonon dispersion and scattering phase space have already been accounted for in the first-principles calculation,³ it raises a question whether the measurement results for thin graphite can be attributed to such confinement effects. In addition, it is important to clarify the thickness dependence that can be produced by phonon Poiseuille flow in mesoscale graphite samples with a thickness more than a few nanometers and less than 1 mm.

The large variation of the thermal conductivity values and trends measured for thin graphite and FLG samples reflect not only the different sample qualities but also the challenges in thermal transport measurements of these small-size, high-thermal conductivity samples. Micro-Raman thermometry measurements of FLG are subject to limited temperature sensitivity, large uncertainty in the optical absorption, and artifacts caused by strain and local non-equilibrium among different phonon polarizations inside the tightly focused laser spot.^23,24 Meanwhile, contact thermal resistance has yielded uncertainties in the results obtained by two-probe thermal transport measurements with electro-thermal microbridge devices, where the complicated microfabrication processes can contaminate the sample surfaces.²⁵ Similar contact resistance error exists in thermocouple-based steady-state measurements of high-thermal conductivity materials, especially mesoscale samples, where the heat loss through the thermocouples can be appreciable compared to heat conduction in the small-cross section samples.²⁶ 

This article attempts to clarify the thickness-dependent thermal conductivity of thin graphite via both numerical solutions of the Peierls–Boltzmann equation (PBE) in three-dimensional (3D) real space and four-probe thermal transport measurements. By assuming the phonon dispersion of graphite and accounting for both phonon–phonon and phonon–surface scattering processes in the anisotropic structure, a deviational Monte Carlo (MC) calculation shows decreased thermal conductivity with decreased thickness from 65 μm to several nanometers, as expected. The calculated temperature dependence and the heat flux profile are used to identify the thickness range for phonon Poiseuille flow to occur in thin graphite. The thickness dependence obtained by the MC simulation is opposite to the recent thermocouple measurement results of thin graphite samples thicker than 8.5 μm.²² However, the MC results with diffuse surface scattering and graphite phonon dispersion are rather close to four-probe thermal transport measurement results of ultrathin graphite samples (UTG) in the thickness range between 13 and 23 layers, which are prepared in this work with an improved method to reduce surface contamination.

In the prior first principles based theoretical calculation of the thermal conductivity of defect-free multi-layer graphene (MLG),³ the unit cell contains the entire MLG thickness to yield a strictly 2D reciprocal space, while a 3D reciprocal space with a small dimension along the cross-plane direction is a suitable representation only when the unit cell is repeated along the thickness direction to yield a periodic superlattice structure. Interlayer interactions in the MLG unit cell result in additional optical phonon polarizations in the 2D Brillouin zone compared to monolayer graphene. These optical polarizations represent essentially non-propagating standing waves along the thickness direction. Such standing waves appear as a result of interference between a lattice wave and its reflection by the specular surface of MLG thinner than the phonon coherent length. As the thickness exceeds the coherence length,²⁷ modes propagating along the thickness direction emerge and replace the standing modes, giving rise to the graphite phonon dispersion at a sufficient thickness. The coherence length is expected to be larger than the phonon wavelength and smaller than the scattering mean free path that governs the lattice thermal conductivity. As the phonon mean free path appears to be about 20 nm on average¹⁹ and may reach about 200 nm^28,29 for heat-carrying phonon modes propagating along the c axis of graphite, the bulk graphite phonon dispersion is expected to be applicable and surface defects would be required to yield a thickness-dependent thermal conductivity in thin graphite samples with a micrometer scale thickness. For UTG with a thickness smaller than 20 nm, a better understanding of the phonon coherence length is required to evaluate the size effects on the phonon dispersion. For different superlattice structures based on either argon or silicon with very different bulk phonon MFPs,³⁰ the coherence length is calculated to be similar and just a few nanometers, somewhat larger than the wavelength of heat-carrying phonons and much smaller than the MFPs of heat-carrying phonons in silicon.³¹ Moreover, the calculated basal-plane thermal conductivity for defect-free MLG approaches the graphite value when the MLG thickness exceeds just five atomic layers.³ Based on this result, treating the interlayer interaction as either propagating modes along the c axis in a 3D graphite phonon dispersion or as additional optical modes in a 2D Brillouin zone do not yield detectable difference in the calculated in-plane thermal conductivity for MLG thicker than five layers.

Theoretical calculations of phonon transport to examine the size effect in three-dimensional finite-size samples with surface defects require the solution of the PBE in real space. Past studies have solved the PBE in real space for steady-state thermal conductivity, but either in low-dimensional materials^32–34 or under relaxation time approximation of the internal phonon scattering processes.^35–37 The exact solution with a full scattering matrix in three-dimensional reciprocal space and three-dimensional real space was not obtained until a very recent work that solved the PBE in all seven dimensions including time domain.¹⁴ 

To investigate the impact of sample dimension and all boundary scatterings of graphite, we solved the PBE in the three-dimensional thin graphite with inputs from first-principles calculations. The steady-state PBE with a full scattering matrix is

where i and j indicate a phonon state, f is the phonon distribution function, $v$ is the phonon group velocity, $r$ is the position vector in the three-dimensional real space, and $Cij$ is the element of scattering matrix $C$⁠. The details of first-principles calculations of phonon dispersion and scattering matrix can be found in a recent literature.¹⁴ The PBE was solved by a deviational MC method, which has been used to study steady-state thermal transport in two-dimensional graphene ribbons^32–35 and transient wave-like thermal transport in three-dimensional graphite.¹⁴ The details of the algorithm can be found in these literatures. Here, we implemented the MC method for thermal conductivity calculations of three-dimensional crystalline graphite.

In the MC calculation, hot and cold reservoirs were attached to the left and right surfaces in the length direction, respectively. Fully diffuse boundary scattering was applied for the side edges of a sample with a finite width (W). For the top and bottom surfaces of thin graphite with a finite thickness (t), the specularity parameter (p) vales were chosen in the range between 0 for fully diffuse scattering and 1 for fully specular reflection.

As a comparison, a solution of the PBE was obtained with the use of an iterative scheme based on the assumption of homogenous f distribution on each cross section to approximate $∇rfi$ as $(dfi0/dT)∇rT$⁠, where $fi0$ is the Bose–Einstein distribution. In this calculation, the boundary scattering rates were taken as the phenomenological models for the diffusive (D) and ballistic (B) phonon transport limits, respectively, according to^38–40 

We start by examining the thickness dependence of the thermal conductivity of thin graphite samples with infinite length and width, and thickness of 8.5 and 65 μm, respectively, which are the two smallest values reported in the recent thermocouple measurement.²² With the phonon dispersion taken as the graphite dispersion, discretization of the rather thick samples into a large number of grids is necessary to obtain sufficient fine grids for the MC result to converge to the expected graphite value when all the surfaces are assumed to be specular. For diffuse surfaces, the thermal conductivity calculated with the graphite dispersion and 25 × 25 × 5 grids in reciprocal space decreases with reduced thickness at temperatures below 150 K, as shown in Fig. 1(a). Above 150 K, the thermal conductivity converges to bulk values calculated for an infinitely large sample and measured for HOPG.⁴² The calculated thickness dependence is negligible, suggesting negligible surface scattering effect in samples thicker than 8.5 μm. In Fig. 1(b), we show the thermal conductivity with thicknesses from a few nanometers to several hundreds of nanometers.

A decreasing thermal conductivity with decreasing thickness is observed at all temperatures, which is shown clearly in Fig. 1(c) and opposite to the trend reported by the recent thermocouple measurement.²² Even with diffuse top and bottom surfaces, in addition, the calculation has yielded much higher thermal conductivity peak values and lower peak temperatures than the thermocouple measurement results. The presence of grain boundaries in the measurement samples are expected to increase the peak temperature and suppress the peak magnitude without increasing the thermal conductivity above that for a single crystal sample of the same thickness, opposite to the behavior shown for the 8.5 μm thick sample in Fig. 1(a). We would like to point out that the higher thermal conductivity for thinner samples reported in Ref. 22 is yet to be properly supported by established theory. The paper connects their observation to the larger size of the first Brillouin zone along the c axis of the thinner samples and, thus, weaker U-scattering involving the reciprocal lattice vector along the c axis. However, the bulk dispersion should still be applicable to graphite samples thicker than a few micrometers. Moreover, reduced or increased U-scattering along the c axis does not affect the thermal transport in the a–b plane, which was discussed in the recent study,^12,41 as it does not destroy the phonon momentum along the heat flow direction.

To investigate whether phonon Poiseuille flow can occur in this thickness range in an intermediate temperature range, we calculate the temperature dependence of thermal conductivity. Figure 2 shows the temperature exponent as a = d(lnκ)/d(lnT) for the three samples at different temperatures compared to the ballistic limit. At 50 K, a calculated for the 8.5 μm sample is much larger than the ballistic limit, indicating the thermal conductivity increases with temperature faster than that in the ballistic case. This behavior is a typical characteristic of phonon Poiseuille flow and agrees with the heat flux profile calculated in Fig. 3 where we show the heat flux across the thickness for the 8.5 μm thick graphite sample at 50 and 100 K. A clear parabolic distribution is observed at 50 K indicating the existence of phonon Poiseuille flow that suppresses the lattice thermal conductivity. As a comparison, the heat flux is uniform at the center at 100 K due to the increased phonon–phonon scattering compared to boundary scattering. To better understand the size dependence of phonon Poiseuille flow, we show the phonon MFPs and the projected MFPs on the in-plane and cross-plane directions in Fig. 4. Due to the high anisotropy in graphite, the z components of MFPs are the relevant length scales for comparison with the sample thickness.⁴³ Both the overall MFPs and the z components by normal scattering are around several micrometers for heat-carrying phonons, in the same range of the thickness of 8.5 μm sample but much smaller than 65 μm.

We further extend the calculation to a 23-layer UTG sample with the length (L), width (W), and thickness (t) being 5, 1.75, and 7.7 nm. The UTG sample was discretized into 50 × 10 × 5 grids in the real-space domain. The first Brillouin zone of reciprocal space was discretized with a 15 × 15 × 3 grid. The time domain was discretized with time steps from 7 to 1 ps at temperatures from 80 to 400 K.

Figure 5(a) includes the calculated MC results for the UTG sample with bulk phonon dispersion, diffuse side edges, and different specularity parameters for the top and bottom surfaces. For fully diffuse scattering at the top and bottom surfaces (p = 0), the MC calculation results are well below those calculated for the bulk samples and the thin graphite samples with 8.5 and 65 μm thickness. As p increases toward 1 for specular scattering at the top and bottom surfaces, the calculated thermal conductivity increases considerably. As a comparison, for samples with infinite length and diffuse scattering at the top and bottom surfaces, expanding the finite side edges to infinite width only increases the thermal conductivity slightly. Moreover, a graphite sample with the same 5 μm length but infinitely large cross section shows only a small increase compared to the UTG sample with specular top and bottom surfaces but diffuse side edges (p = 1), and a small decrease compared to the infinitely large bulk sample. These comparisons reveal that the scattering at the side edges and the two ends yield a much smaller reduction in thermal conductivity compared to diffuse scattering at the top and bottom surfaces in the UTG.

Figure 5(b) shows that the MC calculation results fall within the two iterative solutions calculated with the phenomenological boundary scattering models for the ballistic and diffusive phonon transport limits in the UTG. As shown in Fig. 6, the exponent of the temperature dependence of the calculated thermal conductivity of the UTG sample with either diffuse or specular top and bottom surfaces is still smaller than that of the ballistic limit. Here, phonon Poiseuille flow does not occur in the UTG sample because the normal scattering MFPs are larger than both the 1.75 μm width along the in-plane direction and the 7.7 nm thickness for the cross-plane direction. The optimal width and thickness for phonon Poiseuille flow in graphene and graphite are expected to be several micrometers at temperatures below 100 K, based on the thickness dependence results calculated here and several previous studies of the width dependence.^{33,34,37,41,44}

Due to the aforementioned challenges associated with contact resistance errors and heat loss into thermometers in steady-state thermal conductivity measurements of high-thermal conductivity mesoscale samples, there exist insufficient thermal conductivity measurement data for comparison with the current and past theories of phonon transport in thin graphite samples. Although the recently reported four-probe thermal transport measurement method has addressed these measurement errors,^45,46 the method is limited to samples with a thermal resistance that is comparable to the thermal resistance of the suspended thin film thermometer lines. Consequently, this method can be used to measure only UTG samples with a much smaller cross section and larger thermal resistance than thin graphite samples with a thickness in the micrometer scale. Therefore, only UTG samples have been measured in this work to provide experimental data for comparison with the calculation results for such UTG samples, leaving measurements of thin graphite samples for a further experimental study to address the discrepancy between the current calculation results and the existing thermocouple measurement results.

The UTG samples were exfoliated from different graphite sources on top of the surface of a 300 nm thick SiO₂ grown on a silicon substrate. The graphite sources include natural graphite (NG) and commercial HOPG from both SPI supplies (HOPG 1) and TipsNano (HOPG 2). The thickness of the exfoliated UTG flakes was determined by atomic force microscopy (AFM) measurements. Figure 7(a) shows a prefabricated SiO₂ beam transferred with the use of a sharp tungsten tip on top of the exfoliated flake. The oxide beam acted as a hard mask during oxygen plasma patterning of the graphite flake sample. Following the patterning, a layer of polymethylmethacrylate (PMMA) was spun onto the substrate to act as a carrier film for subsequent transfer of the patterned flake and oxide beam hard mask, which prevented the PMMA from contacting the top surface of the graphite flake. For the transfer, the PMMA film along with the sample and oxide beam stack was released from the silicon substrate via dilute hydrofluoric acid (HF) etching of the oxide film under the graphite flake. The released PMMA film was then transferred onto multiple Pd/Cr/SiN_x thermometer lines that were patterned on top of another Si substrate, as shown in Fig. 7(b). Electron beam lithography was used to open windows in the PMMA at the intersection of the oxide beam and the underlying metal resistance thermometer lines. Diluted HF was used to etch the exposed SiO₂ beam through the open window in the PMMA, followed by evaporation of Cr and Pd and removal of the PMMA in acetone to lift off Pd/Cr anchors to fix the graphite flake onto the thermometer lines. The thermometer lines and the graphite sample were then suspended by etching of the underlying silicon in tetramethylammonium hydroxide (TMAH) or potassium hydroxide (KOH). The oxide beam was subsequently removed in HF. The measurement device was dried in a critical point dryer. Figure 7(c) shows the suspended UTG sample on the four-probe microdevice. The suspended sample was loaded in the evacuated sample space of a cryostat for multi-probe thermal transport measurements. The measurement was performed by resistively heating each single metal thermometer line while simultaneously measuring the electrical resistances of all thermometer lines, as described in a previous report.⁴⁵ 

Upon the completion of the thermal measurements, Raman spectroscopy, scanning electron microscopy (SEM), and atomic force microscopy (AFM) were used to characterize the UTG samples directly on the suspended devices. The SEM measurements obtained the lateral dimension of each suspended segment, as summarized in Table I. The Raman excitation laser wavelength was 532 nm. Subsequently, a tungsten (W) tip attached to a micromanipulator under an optical microscope was used to transfer one suspended segment of sample NG 1-b onto a gold-coated Si substrate for scanning tunneling microscopy (STM) characterizations of the sample surface. The STM measurement was carried out inside the ultrahigh vacuum (UHV) chamber of an Omicron scanning tunneling microscope. The vacuum level was maintained below 9.0 × 10⁻¹⁰ Torr. During the STM measurement, a sharp W tip prepared by electrochemical etching was scanned at a constant current of 50 pA under 0.1 V voltage applied between the tip and the sample. All surface characterizations were carried out on the top surface of the suspended UTG sample.

The measured thermal conductivity of UTG samples peaks at a temperature near 275 K with the value ranging from 500 to 706 W m⁻¹ K⁻¹, regardless of sample sources, in Fig. 5(b). The obtained thermal conductivity values at room temperature are close to a quarter of the basal-plane value reported for high-quality pyrolytic graphite (PG), for which the peaks appear at a lower temperature near 100 K than that observed for the UTG samples.^19,42 Except for NG2 that shows the lowest thermal conductivity in the six samples measured here, the thermal conductivity peak measured for the other UTG samples is slightly increased and shifted to a lower temperature compared to the reported 480 ± 21 W m⁻¹ K⁻¹ peak value at near 300 K for a seven-layer FLG sample that was transferred onto a four-probe device without the use of the SiO₂ dry mask to separate the PMMA and the FLG sample.⁴⁵ In comparison, a peak amplitude of about 1000 W m⁻¹ K⁻¹ appears at near 230 K for a reported 27-layer thick sample supported on SiO₂.¹⁹ Both the peak magnitude and position reveal the dominance of extrinsic phonon scattering processes compared to intrinsic phonon–phonon scattering in the UTG samples.

The UTG measurement results do not reveal a clear thickness dependence. However, both the magnitudes and temperature dependence of the UTG samples are similar to the MC calculation results for a 23-layer sample with the bulk phonon dispersion and diffuse top and bottom surfaces. As the MC calculations reveal that phonon scattering at the side edges and two ends are insufficient to cause the observed pronounced effect, defects on the surfaces or in the interior of the UTG samples need to be examined.

The Raman spectra (Fig. 8) measured on the suspended UTG segment does not show clear D peak associated with defect-mediated Raman scattering. In comparison, AFM measurements performed directly on the top surface of the suspended segment reveals the presence of nanoparticles with a height over 3 nm on the top surface of NG 1-b, as shown in Fig. 9(a). The areal coverage of nanoparticles is about 6.5% based on the AFM images. These nanoparticles are likely placed on both top and bottom surfaces during the sample suspension step, the only step where the top surface had not been protected by the SiO₂ beam. The three-dimensional height image of the selected area [Fig. 9(b)] showed a root-mean-square height fluctuation of 464.4 pm without any evidence of kinks or steps. This height fluctuation is likely due to the vibration of the suspended segment during the AFM measurement.

In addition, the STM measurements of the UTG sample transferred to the Au-coated Si substrate can help to examine whether there exists a polymeric residue or atomic disorders on the sample surface. Figures 10(a) and 10(b) show the STM height and current response over a 400 × 400 nm² scan area on the UTG sample. Some long-range modulation about 100 nm with height variation less than 2 nm was observed in the STM images of the supported UTG. The modulation was likely caused due to the surface feature of Au-coated Si substrate since the AFM measurement results on the same original suspended UTG segment do not show similar modulation. Figure 10(c) shows the STM height image for a 2 × 2 nm² scan area. Although the spatial and height resolutions of the room-temperature STM measurements were limited by noise sources including vibration and thermal drift, most of the scanned area do not exhibit clear evidence of atomic disorder or point defect. Nevertheless, the scanned area shows several dark or bright spots, which are potentially caused by point defects based on the prior STM studies on point defects on graphite surfaces.^47,48 The bottom surface of the UTG sample was not accessible for the STM and AFM measurements. Because the bottom surface was exposed to an acetone solution after the PMMA was dissolved in the solution, the bottom surface was likely contaminated by polymer residue. In addition to point defects on the surfaces, contamination of the sample surface by nanoparticles and polymer residue potentially results in partially diffuse surface scattering at the top and bottom surfaces, respectively. The lack of clear thickness dependence is attributed to a variation of the surface specularity parameter as well as interior defects that have not been clearly revealed by the Raman and STM measurement results.

The presented deviational MC analysis has calculated the spatial dependent phonon heat flux in thin graphite as thick as 65 μm and as thin as 3.35 nm for a 10-layer graphene. Based on the prior first-principles calculation that shows a saturation to the graphite thermal conductivity for few-layer graphene samples thicker than about five layers, bulk graphite phonon dispersion has been used in the MC calculation. For thin graphite samples with a thickness of 8.5 μm and above, the use of bulk graphite phonon dispersion should be accurate because the phonon coherent length is not expected to reach this rather large thickness to yield non-propagating standing waves along the cross-plane direction, whereas a reduced Brillouin zone along the cross-plane direction may be achieved only in superlattice structures instead of a homogeneous thin film. With the bulk graphite phonon dispersion, thickness-dependent thermal conductivity can only be caused by diffuse scattering by surface defects, which decreases the basal-plane thermal conductivity with decreased thickness and produces an opposite trend as the recent thermocouple measurements.²² Based on the MC calculation, phonon Poiseuille flow occurs in finite-size graphite near 50 K when the thickness is around several micrometers and comparable to the normal scattering MFPs along the thickness direction. Remarkably, MC calculation with the bulk phonon dispersion and diffuse surface scattering yields results close to the measured thermal conductivity for UTG samples with a thickness close to 7.7 nm or 23 layers, where phonon scattering at the side edges of the 1.75 μm wide sample and two ends of the 5 μm long UTG are insufficient to suppress the thermal conductivity to the measured value. While the observed contamination and potentially point defects on the surfaces may provide mechanisms for partially diffuse surface scattering, whether strain or a modified phonon dispersion contributes to the suppressed thermal conductivity in the suspended UTG samples provides an intriguing question for further theoretical studies. Instead of providing definite answers to these outstanding questions, this reexamination serves to highlight the need of further experimental and theoretical studies to address the discrepancy between existing theoretical and experimental results in thin graphite and to minimize extrinsic defect-mediated processes to unmask the intrinsic low-dimensional phonon transport behaviors in ultrathin graphite and few-layer graphene samples.

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

The collaboration between the theory and experiments was supported by National Science Foundation (NSF) (Award Nos. CBET-1705756 and CBET-1707080). X.L. acknowledges support for theoretical calculations from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. H.L. was supported by U.S. Office of Naval Research MURI (Grant No. N00014-16-1-2436). The calculations used resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. The experimental work was partly done at the Texas Nanofabrication Facility supported by NSF grant NNCI-2025227.

The authors have no conflicts of interest to disclose.

X.L., H.L. and E.O. contributed equally to this work.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.