Controlling Thermal Conductivity of Two-dimensional Materials via Externally Induced Phonon-Electron Interaction

TL;DR

This paper investigates how external electrostatic gating can modulate phonon-electron interactions in 2D materials, significantly reducing their thermal conductivity and enabling potential thermal switching applications.

Contribution

It introduces an ab initio approach and an analytical model to understand and control phonon-electron scattering in 2D materials, a previously unexplored area.

Findings

Over 40% reduction in thermal conductivity in silicene with charge carrier concentration of 10^13 cm^-2

Demonstrates feasibility of thermal conductivity control via electrostatic gating in 2D materials

Provides guidelines for designing 2D materials with enhanced phonon-electron scattering

Abstract

Phonon scattering by electrons, or "phonon-electron scattering", has been recognized as a significant scattering channel for phonons in materials with high electron concentration, such as thermoelectrics and nanoelectronics, even at room temperature. Despite the abundant previous studies of phonon-electron scattering in different types of three-dimensional (3D) bulk materials, its impact on the phonon transport, and thus the heat transfer properties, of two-dimensional (2D) materials has not been understood. In this work, we apply ab initio methods to calculate the phonon-electron scattering rates in two representative 2D materials, silicene and phosphorene, and examine the potential of controlling the thermal conductivity of these materials via externally induced phonon-electron scattering by electrostatic gating. We also develop an analytical model to explain the impact of reduced dimensionality and distinct electron and phonon dispersions in 2D on phonon-electron scattering processes. We find that over 40\% reduction of the lattice thermal conductivity can be achieved in silicene with an induced charge carrier concentration in the range of $10^{13}~cm^{-2}$, which is experimentally achievable. Our study not only generates new fundamental insights into phonon transport in 2D materials but also provides practical guidelines to search for 2D materials with strong phonon-electron scattering for potential thermal switching applications.

Figures

Equations

$$\frac{1}{{\tau }_{\mathbf{q}\nu }^{\mathrm{ep}}}=-\frac{2\pi }{\hslash }\sum _{mn,\mathbf{k}}|g_{mn}^{\nu }(\mathbf{k},\mathbf{q}){|}^2(f_{n\mathbf{k}}-f_{m(\mathbf{k}+\mathbf{q}+\mathbf{G})})\delta ({ϵ}_{n\mathbf{k}}-{ϵ}_{m(\mathbf{k}+\mathbf{q}+\mathbf{G})}-{\omega }_{\mathbf{q}\nu })$$

$$\frac{1}{{\tau }_{\mathbf{q}\nu }^{\mathrm{ep}}}=\frac{D^2(2\pi m^{\ast }{)}^{\frac{1}{2}}}{\rho {\alpha }^{\frac{1}{2}}(k_{\mathrm{B}}T{)}^{\frac{3}{2}}}{e}^{-\frac{{\hslash }^2(\frac{1}{2}+\frac{\alpha m^{\ast }}{\hslash }{)}^2{\omega }_{\mathbf{q}\nu }}{2m^{\ast }\alpha k_{\mathrm{B}}T}}n(E_{\mathrm{F}}){\omega }_{\mathbf{q}\nu }^{\frac{1}{2}},$$

$$\frac{1}{{\tau }_{\mathbf{q}\nu }^{\mathrm{ep}}}=\frac{D^2(2\pi m^{\ast }{)}^{\frac{1}{2}}}{\rho v_{\mathrm{s}}(k_{\mathrm{B}}T{)}^{\frac{3}{2}}}{e}^{-\frac{{\hslash }^2{\omega }_{\mathbf{q}\nu }^2}{8m^{\ast }{v}_{\mathrm{s}}^2{k}_{\mathrm{B}}T}-\frac{\hslash {\omega }_{\mathbf{q}\nu }}{2k_{\mathrm{B}}T}-\frac{m^{\ast }{v}_{\mathrm{s}}^2}{2k_{\mathrm{B}}T}}n(E_{\mathrm{F}}){\omega }_{\mathbf{q}\nu },$$

$$\mathrm{TDOS}(E_{\mathrm{F}})={\int }_{-\infty }^{+\infty }\mathrm{DOS}(E)f^{\prime }(E,E_{\mathrm{F}})dE,$$

Full text

Electron-phonon interactions play a major role in determining the electronic properties of materials since they are the major contributors to electrical resistance and also mediate conventional superconductivity1, 2. Due to these reasons, the influence of electron-phonon interaction on the transport of electrons has been intensively studied and well understood. On the other hand, the scattering of phonons due to electron-phonon interactions (hereafter “phonon-electron scattering”) and its impact on thermal transport of solids have received limited research interest, due to the long belief that it is only important at cryogenic temperatures3, 4, 5. The main reason is that most of the previous studies and practical interests were limited to devices with a low or moderate electron concentration, typically below $10^{19}\ \mathrm{cm}^{-3}$. Recent technological developments have led to important applications with electron concentration as high as $10^{20}$ to $10^{21}\ \mathrm{cm}^{-3}$, such as in heavily-doped thermoelectric materials6 and nanoelectronic devices7. In this regime, however, the impact of phonon-electron scattering on thermal transport in largely unknown. Recently, Liao et al. used ab initio calculations8 to show that the lattice thermal conductivity of silicon with a high electron concentration can be suppressed by as much as 50% even at room temperature due to phonon-electron scattering9. Significant suppression of phonon propagation by phonon-electron scattering in silicon at room temperature was subsequently verified experimentally using ultrafast photoacoustic spectroscopy10. Moreover, ab initio calculations of phonon-electron scattering have also been carried out in bulk metals11, 12, providing new insights into the details of coupled transport of phonons and electrons.

With the rapid advancement of nanotechnology, 2D materials have become star candidates for a wide range of applications, e.g. transistors, optoelectronics and energy harvesting devices. Currently, theoretical and experimental understanding of electron-phonon interaction in 2D materials has been limited to its influence on electrons and charge carrier mobility13, 14, 15. The effect of phonon-electron scattering on the phonon frequency renormalization of 2D materials has been reported16, but its effect on phonon lifetime and transport properties is still lacking. Given the paramount importance of phonon transport and thermal management for device performance in these applications, it is desirable to gain an in-depth understanding of phonon-electron scattering and its impact on phonon transport in 2D materials. From a fundamental point of view, phonon-electron scattering in 2D materials is expected to be qualitatively different from that in 3D bulk materials, due to factors including reduced dimensionality and thus altered scattering phase space, dominant normal phonon scattering and hydrodynamic phonon transport17, 18, new symmetries and the associated scattering selection rules19, distinct electron and phonon dispersion relations (e.g. massless Dirac fermions and quadratic flexural phonons) and different dielectric screening behavior for polar materials20. Therefore, it is of fundamental interest to explore the behavior of phonon-electron scattering in 2D materials. From a practical point of view, thanks to the possibility of inducing a high concentration of electrons or holes in 2D materials via electrostatic gating15, 21, efficient phonon-electron scattering in 2D materials can potentially enable the development of fast thermal switches whose thermal conductivity can be tuned by applying an external electric field. This concept is illustrated in Fig. 1a. Reversible and wide-range control of the thermal conductivity of solids using external fields is highly desirable in diverse fields22. Previously, a variety of approaches have been explored experimentally to reversibly control the thermal conductivity of solids. For example, the metal-insulator phase transition23 and the associated crystal lattice change, e.g. of vanadium dioxide, has been proposed as a potential mechanism of thermal switching, but only a small contrast was achieved experimentally24. In another example, a thermal conductivity tuning of 11% was demonstrated by modifying the ferroelectric domain structure using an external electric field25. Reversible electrochemical intercalation of ions has been demonstrated to change the thermal conductivity of layered materials26, but requires operation in a liquid electrolyte. Here we provide another possible scheme. With the development of solid-state electrolytes as the gate dielectric, inducing a charge concentration as high as $10^{14}\ \mathrm{cm}^{-2}$ in 2D materials15, 21 has recently been demonstrated. In principle, this high density of charge carriers can efficiently scatter phonons and largely reduce the lattice thermal conductivity while an optimum charge carrier density can be selected to balance the increase of the electronic thermal conductivity.

In this work, we use ab initio electron-phonon interaction calculations8 to examine the practicality of this mechanism as a means to realize efficient thermal switching of 2D materials. We choose silicene and phosphorene as model systems due to their wrinkled crystal structures that break the out-of-plane mirror symmetry, as shown in Fig. 1b. It has been theoretically shown that the out-of-plane mirror symmetry, e.g. of graphene, prohibits the first-order interaction of electrons with the out-of-plane flexural phonons27, 28, which are the major heat carriers in 2D materials19. Therefore, 2D materials without the out-of-plane mirror symmetry should be better candidates to possess strong phonon-electron scattering. Another reason to choose silicene and phosphorene is their distinct electronic structures: silicene is a Dirac semimetal with massless electrons and holes near the intrinsic Fermi level29 and phosphorene is a semiconductor with a sizable band gap and highly anisotropic carrier effective masses30, 14. Therefore, a comparative study will help elucidate the impact of electronic structure on phonon-electron scattering and provide guidelines in future search for 2D materials with desirable phonon-electron scattering properties. We find that the lattice thermal conductivity of silicene and phosphorene depends strongly on the induced charge carrier concentration, as illustrated in Fig. 1c. In particular, in p-type silicene, over 40% reduction of the lattice thermal conductivity can be achieved by a charge carrier density of $10^{13}\mathrm{cm}^{-3}$, whereas in phosphorene, 10% change of the lattice thermal conductivity can be expected at a similar charge carrier density. We provide in-depth analysis of these results in the following sections.

The scattering rate of a phonon mode with wave vector $\mathbf{q}$ and branch label $\nu$ due to phonon-electron scattering is given by the Fermi’s golden rule in the form2, 9:

[TABLE]

where $\tau^{\mathrm{ep}}_{\mathbf{q}\nu}$ is the corresponding phonon relaxation time, $g^{\nu}_{mn}(\mathbf{k},\mathbf{q})$ is the scattering matrix element connecting the initial electronic state with band index $n$ and wave vector $\mathbf{k}$ with the final electronic state with band index $m$ and wave vector $\mathbf{k}+\mathbf{q}+\mathbf{G}$ ($\mathbf{G}$ is a reciprocal lattice vector) due to the conservation of crystal momentum, $f_{n\mathbf{k}}$ is the Fermi-Dirac distribution, $\epsilon_{n\mathbf{k}}$ is the electron energy, and $\omega_{\mathbf{q}\nu}$ is the phonon frequency. The $\delta$ function imposes energy conservation during the scattering process. From Eq. 1, there are two major factors that contribute to the phonon-electron scattering rate: the scattering matrix elements that reflect the strength of the coupling between the electronic and lattice degrees of freedom and the energy-momentum conservation conditions that determine the number of potential scattering channels. To carry out the calculation in Eq. 1, the electronic band structure, the phonon dispersion and electron-phonon scattering matrix element are required. We use the density functional theory (DFT) based method to calculate these ingredients. The electronic band structure is calculated using the Quantum ESPRESSO package31. A mesh grid of $30\times 30\times 1$ in the first Brillouin zone is adopted for both materials with norm-conserving pseudo-potentials. The kinetic energy cutoff for wavefunctions is set to 40 Ry. The kinetic energy cutoff for charge density and potential is set to 160 Ry. The total electron energy convergence threshold for self-consistency is $\rm 1\times 10^{-10}~{}Ry$. Applying density functional perturbation theory (DFPT) implemented in the same package, the phonon dispersion and the electron-phonon matrix elements are calculated on a coarse mesh of $6\times 6\times 1$ and $4\times 4\times 1$ for silicence and phosphorene, respectively. The results are subsequently interpolated onto a fine mesh of $60\times 60\times 1$ for silicene and $64\times 64\times 1$ for phosphorene using a maximally localized Wannier function based scheme as implemented in the EPW package32. The calculated electronic band structure and phonon dispersion agree well with the literature and are given in the Supporting Information.

Figure 2 shows the scattering rates of different phonon modes due to phonon-electron scattering in phosphorene. Both n-type and p-type phosphorene with different charge carrier concentrations are investigated. Several observations can be made compared to the results in 3D materials, such as bulk silicon9. While the phonon scattering rates scale linearly with the carrier concentration in both silicon and phosphorene, a major difference is the phonon frequency dependence of the scattering rates for the low-frequency acoustic phonons. In 3D semiconductors like silicon, the linear phonon dispersion and parabolic electronic bands give rise to the linear dependence of the phonon scattering rates on the phonon frequency for the low-frequency acoustic phonons9, whereas in phosphorene, the frequency dependence is more complicated: the scattering rates first rise with phonon frequency below 0.3 THz and then decrease. The decreasing trend of the phonon scattering rates in phosphorene indicates that the major heat-carrying phonons, typically around 1 THz, are less effectively scattered by electrons than phonons with longer wavelengths. To qualitatively understand the distinct phonon frequency dependence of the phonon-electron scattering rates in phosphorene, we develop an analytical model for the phonon-electron scattering rates of low-frequency acoustic phonons with the deformation potential approximation, where the electron-phonon matrix elements $|g^{\nu}_{mn}(\mathbf{k},\mathbf{q})|^{2}$ are replaced by $\frac{\hbar D^{2}q^{2}}{2m_{0}\omega_{\bf{q}\nu}}$ with $D$ being the constant deformation potential33 and $m_{0}$ the atomic mass. The deformation potential is the change of the electron energy near the band edges in response to a static strain of the lattice, corresponding to the coupling strength of electrons with long-wavelength acoustic phonons. Using the phonon dispersion relations in 2D, including the quadratic dispersion for out-of-plane flexural phonons, and the 2D electronic band structure, we derive the following analytical equations for the phonon-scattering rates of ZA phonons in 2D semiconductors in the nondegenrate regime:

[TABLE]

where $m^{*}$ is the average electron or hole effective mass, $\rho$ is the areal mass density, $\alpha$ is the coefficient for the quadratic ZA phonon dispersion ($\omega=\alpha q^{2}$), $k_{\mathrm{B}}$ is the Boltzmann constant and $n(E_{\mathrm{F}})$ is the carrier concentration corresponding to a certain Fermi level $E_{\mathrm{F}}$. The detailed derivation is given in the Supporting Information. Similarly, we also derive an analytical expression for the phonon scattering rates of TA and LA phonon modes with linear dispersion relations given by:

[TABLE]

where $v_{\mathrm{s}}$ is the speed of sound for the TA and LA phonon modes. These analytical results are plotted in Fig. 2 as solid lines and agree very well with the ab initio calculations, and confirm that the different phonon frequency dependency originates from the difference in dimensionality and the electron and phonon dispersion relations. These formulas can be used to estimate the level of phonon-electron scattering rates for long-wavelength acoustic phonons in 2D semiconductors. In Fig. 2(d) and (h), we also report the phonon-mode-resolved phonon-electron scattering rates in n-type and p-type phosphorene, where each data point in the plots represents a phonon mode with given frequency and wavevector magnitude and the color of the data point encodes the phonon-electron scattering rate of this mode. In both n-type and p-type phosphorene, phonons near the zone center are much more strongly scattered by electrons. This behavior is also observed in bulk silicon9 and is expected due to the small energy scale of phonons compared to that of electrons, so that the small-wavevector phonons have a larger chance of meeting the energy-momentum conservation conditions during the scattering process. Furthermore, some of the zone boundary phonons, particularly in p-type phosphorene, are also significantly scattered by electrons. This is related to the strong anisotropy and large carrier effective mass along the zigzag direction in phosphorene30, 14, which leads to large radius of the electron and hole pockets near the band edges such that the electrons can be scattered across the carrier pockets by the zone boundary phonons with a large momentum change but a small energy change.

In Fig. 3, we show the calculated scattering rates of phonons due to phonon-electron scattering in silicene. The main difference between silicene and phosphorene is their electronic structures: phosphorene is a wide-gap semiconductor and silicene is a Dirac semimetal with linear electron bands. The phonon frequency dependence of the scattering rates of ZA mode shows similar behavior as that in phosphorene, but the scattering rates are higher in silicene at a given carrier concentration in general. As the carrier concentration increases, the phonon scattering rates also exhibit a trend of saturation as the Fermi level moves deeper into the bands. Due to the linear electron bands and the degenerate nature requiring the use of Fermi-Dirac distribution functions in Eq. 1, there is no simple analytical expression for the phonon scattering rates as those given for phosphorene (Eq. 2 and Eq. 3) even with the deformation potential approximation. Instead, we integrate Eq. 1 numerically with the deformation potential approximation, and the results are shown as solid lines in Fig. 3, where good agreement with the ab initio results is observed and the saturation of the scattering rates with increasing carrier concentration is captured. Figure 3(d) and (h) display the mode-resolved phonon-electron scattering rates in n-type and p-type silicene. Similarly as in phosphorene, phonons near the zone center are strongly scattered by electrons. Phonons in a section of the Brillouin zone near the boundary are also strongly scattered, which is caused by intervalley scattering of electrons and holes between the two Dirac cones of silicene.

From Figs. 2 and 3, it is clear that the phonon scattering rates due to phonon-electron scattering are significantly higher in silicene than those in phosphorene given similar carrier concentrations. To understand this difference, we analyze the strength of phonon-electron scattering by considering the two major factors that determine the phonon-electron scattering rates as reflected in Eq. 1: the number of potential phonon-electron scattering channels imposed by the energy-momentum conservation conditions and the magnitude of the electron-phonon scattering matrix elements. Given the small energy scale of phonons compared to that of electrons, phonon-electron scatterings are approximately “on-shell” processes, where the initial and the final electron states have similar energy. This observation implies that the number of potential phonon-electron scattering channels is directly related to the electron density of states within the Fermi window, where scatterings can happen between partially occupied energy levels. A quantitative measure is the so-called “thermal density of states” (TDOS)2, which is defined as a function of the Fermi level:

[TABLE]

where $\mathrm{DOS}(E)$ is the electron density of states, and $f^{\prime}(E,E_{\mathrm{F}})$ is the energy derivative of the Fermi-Dirac distribution, which defines the Fermi window. In Fig. 4(a) and (b), the TDOS for both n-type and p-type silicene and phosphorene are shown, which clearly signal that there is a larger available scattering phase space in silicene than phosphorene at the same carrier concentration due to silicene’s semimetallic band structure. In Fig. 4(c) and (d), we further calculate and compare the scattering matrix elements for acoustic phonon modes in silicene and phosphorene given the same initial electron states at $\mathbf{k}=\mathbf{0}$ and at the band extrema. The electron-phonon matrix elements can be calculated as $g^{\nu}_{mn}(\mathbf{k},\mathbf{q})=\Braket{u_{m(\mathbf{k}+\mathbf{q})}}{\Delta_{\mathbf{q}\nu}v^{\mathrm{SCF}}}{u_{n\mathbf{k}}}$, where $u_{m(\mathbf{k}+\mathbf{q})}$ and $u_{n\mathbf{k}}$ are the periodic parts of the Bloch wavefunctions of the corresponding electron states and $\Delta_{\mathbf{q}\nu}v^{\mathrm{SCF}}$ is the perturbation of the electron self-consistent-field potential induced by a phonon mode with wavevector $\mathbf{q}$ and branch index $\nu$8. Physically, the matrix elements measure the sensitivity of the electron potential energy in response to lattice disturbance and can be affected by crystal structure, chemical bonding environment and electron screening effect2. Here we see that the electron-phonon scattering matrix elements in silicene are significantly larger in magnitude than those in phosphorene, signaling that the electron states in silicene are much more sensitive to lattice perturbation and the structure of silicene might thus be less stable than that of phosphorene.

Finally, to compare the contribution of phonon-electron scattering to the lattice thermal conductivity with that of intrinsic phonon-phonon scattering, we calculate the phonon-phonon scattering rates in silicene and phosphorene based on the anharmonic interatomic force constants (IFCs). The anharmonic IFCs are calculated using a frozen-phonon approach with supercell sizes $\rm 5\times 5\times 1$ and $\rm 4\times 4\times 1$ for silicene and phosphorene, respectively. The q-mesh grid $60\times 60\times 1$ is adopted for both silicene and phosphorene. The results are shown in Fig. 5. When the carrier concentration is higher than $1\times 10^{13}\ \mathrm{cm}^{-2}$, the scattering rates of the low-frequency acoustic phonons below a few THz are dominated by phonon-electron scattering, which is expected to have a major impact on the lattice thermal conductivity of these materials. To quantify this effect, we further combine the phonon-phonon and phonon-electron scattering rates using Matthiessen’s rule and solve the phonon Boltzmann transport equation with the ShengBTE package34, 35 to obtain the lattice thermal conductivity as a function of the carrier concentration, which is shown in Fig. 1(c). Due to the ambiguity of defining an effective thickness for 2D materials to calculate their thermal conductivity, here we report the relative reduction of the lattice thermal conductivity due to phonon-electron scattering normalized to the intrinsic values without net induced charge. The intrinsic thermal conductance (product of layer thickness and thermal conductivity) we obtain is $5.3\times 10^{-9}\ \rm{W/K}$ for silicene, $5.8\times 10^{-8}\ \rm{W/K}$ along the zigzag direction of phosphorene and $2.3\times 10^{-8}\ \rm{W/K}$ along the armchair direction of phosphorene, which is in good agreement with literature values36, 37. Our results demonstrate that externally induced phonon-electron scattering can significantly affect thermal transport in silicene and phosphorene, particularly in silicene with high TDOS and strong electron-phonon coupling as reflected in the magnitude of the matrix elements. At a high carrier concentration of $10^{13}\ \rm{cm}^{-2}$, a 4̃0% reduction of the lattice thermal conductivity in p-type silicene can be achieved, indicating the potential use of this mechanism to realize thermal switching devices. One potential concern is the increased electronic thermal conductivity as the induced carrier concentration increases. Due to the lack of available experimental data of electronic thermal conductivity in these materials, we estimate the electronic thermal conductivity using experimentally reported mobility38, 39 and the Wiedemann-Franz law to be $\sim$ 10% and $\sim$ 5% of the intrinsic lattice thermal conductivity in silicene and phosphorene at a carrier concentration of $10^{13}\ \rm{cm}^{-2}$, respectively. It is worth noting that, in addition to inducing charge carriers, the external electric field applied to 2D materials can also cause changes to charge density distribution inside the material that can further reduce the thermal conductivity, as has been studied from first-principles in silicene40. We envision that a combination of these effects will render silicene a promising candidate for thermal switching applications driven by an external electric field.

In summary, we analyze the impact of phonon-electron scattering on phonon transport in two representative 2D materials, silicene and phosphorene, using ab initio calculations. We examine the mode-resolved phonon-electron scattering rates in silicene and phosphorene in detail and compare their behavior to 3D bulk materials. We explain the observed phonon-frequency dependence by developing a semi-analytical model using the deformation potential approximation, based on which we attribute the qualitatively different behavior to distinct phonon and electron dispersion relations in 2D materials as well as the reduced dimensionality. We provide further understanding about significantly different phonon-electron scattering rates in silicene and phosphorene by analyzing the magnitude of the scattering matrix elements and the available scattering phase space quantified by the thermal density of states. We find that the lattice thermal conductivity of silicene can be reduced by over 40% via externally induced phonon-electron scattering, indicating the potential application in thermal switching devices. Our fundamental study also provides guidelines to search for 2D materials with even stronger phonon-electron scatterings for thermal switching applications.

1 Author Contributions

S. Y. and B. L. initiated this project. S. Y. carried out the ab initio calculations. R. Y. developed the semi-analytical models. S. Y. and B. L. wrote the manuscript.

{acknowledgement}

This work is based on research supported by the Army Research Office Young Investigator Program under the award number W911NF-19-1-0060. We acknowledge the use of the Center for Scientific Computing supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC) at University of California, Santa Barbara through NSF awards DMR-1720256 and CNS-1725797. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE) Stampede 2 at the Texas Advanced Computing Center (TACC) through allocation TG-DMR180044. XSEDE is supported by NSF grant number ACI-1548562.