Unusual Enhancement in Intrinsic Thermal Conductivity of Multilayer

Aug 4, 2015 - Due to the significant challenges of k measurements at nano- and microscales, the available data for k of suspended MLG is limited, and ...
0 downloads 11 Views 990KB Size
Subscriber access provided by UNIV OF OTTAWA

Communication

Unusual Enhancement in Intrinsic Thermal Conductivity of Multi-layer Graphene by Tensile Strains Youdi Kuang, Lucas Lindsay, and Baoling Huang Nano Lett., Just Accepted Manuscript • DOI: 10.1021/acs.nanolett.5b02403 • Publication Date (Web): 04 Aug 2015 Downloaded from http://pubs.acs.org on August 6, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Nano Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

Unusual

Enhancement

in

Intrinsic

Thermal

Conductivity of Multi-layer Graphene by Tensile Strains Youdi Kuang,* †‡ Lucas Lindsay,§ and Baoling Huang‡ †

College of Engineering, Shanghai Second Polytechnic University, Pu Dong, Shanghai, China



Department of Mechanical and Aerospace Engineering, The Hong Kong University of Science

and Technology, Clear Water Bay, Kowloon, Hong Kong §

Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge,

Tennessee 37831, USA

ABSTRACT: Using the Boltzmann-Peierls equation for phonon transport approach with the inputs of interatomic force constants from the self-consistent charge density functional tight binding method, we calculate the room-temperature in-plane lattice thermal conductivities k of multi-layer graphene (up to four layers) and graphite under different isotropic tensile strains. The calculated in-plane k of graphite, finite mono-layer graphene and 3-layer graphene agree well with previous experiments. For unstrained graphene systems, both the intrinsic k and the extent of the diffusive transport regime present a drastic dimensional transition in going from monolayer to 2-layer graphene and thereafter a gradual transition to the graphite limit. We find a peak enhancement of intrinsic k for multi-layer graphene and graphite with increasing strain with the

ACS Paragon Plus Environment

1

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 19

largest enhancement amplitude ~40%. Competition between the decreased mode heat capacities and the increased lifetimes of flexural phonons with increasing strain contribute to this k behavior. Similar k behavior is observed for 2-layer hexagonal boron nitride systems. This study provides insights into engineering k of multi-layer graphene and boron nitride by strain and into the nature of thermal transport in quasi-two-dimensional and highly anisotropic systems.

KEYWORDS: Tensile strain, density functional tight binding, thermal conductivity, multi-layer graphene, phonon thermal transport

In contrast to purely two-dimensional (2D) mono-layer graphene, multi-layer graphene (MLG) has a quasi-2D structure with weak coupling between layers.1 Though such coupling decreases the thermal conductivity k compared with mono-layer graphene, experimentally reported k of MLG still range up to 2800 W/m-K,1, 2 far larger than those of most thermally conductive materials. Therefore, MLG is an ideal candidate for thermal management of nanoelectronics 1 and nanocomposites3, and has attracted wide interest in the past decade.1, 2, 4-7 Due to the significant challenges of k measurements at nano- and micro-scales, the available data for k of suspended MLG is limited, and measured values strongly depend on a variety of factors including in-plane size8 and layer number1 of MLG, measurement temperature,5 measurement method2 and sample quality.9 Using an optothermal Raman technique, Ghosh et al.1 characterized the dimensional transition of k in MLG. They observed that k decreases from 2800 W/m-K to 1300 W/m-K at room temperature as the layer number n increases from n=2 to n = 4 but the measured k lie below that of graphite (kgraphite ≈ 2000 W/m-K10) for n = 4 and n ≈ 8. Li’s2 group obtained room-temperature k ≈ 1400 W/m-K and k ≈ 1490 W/m-K for n = 3 for two

ACS Paragon Plus Environment

2

Page 3 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

samples of different width using the direct thermal-bridge measurement technique. Following the same method, Shi’s group9 found that polymer residue from transfer processing has a significant effect, suppressing k of bi-layer graphene. Suppression of k in tri-layer graphene with Au deposition was also observed by Li’s group.2 The k of Bernal stacked and twisted bi-layer graphene11 was reported to have near room temperature values of 1900 W/m-K and 1410 W/mK, respectively. The lower k in Ref. 11 compared to Ref. 1 is attributed to smaller sample size and twisting defects. On the theoretical side, previous molecular dynamics (MD) simulations12 using the Brenner empirical potential13 showed in-plane k decreases with increasing layer number and saturates at n = 4 but with a much lower in-plane k value than that of graphite. Based on the iterative solution of the linearized Boltzmann-Peierls Equation (BPE) with interatomic force constants (IFCs) from the optimized Tersoff potential,14 Lindsay et al.4 found flexural phonons in MLG contribute dominantly to k of finite graphene (less than 10 µm in length) with strong boundary scattering even though a symmetry selection rule in mono-layer graphene is broken; such dominant contribution was also observed by Singh et al.7 using the same potential. However, recent equilibrium MD simulations15 using this potential gave much lower converged room-temperature k (~1000 W/m-K) of mono-layer graphene than experiments16. In contrast to the use of empirical potentials, Paulatto et al.17 explored the anharmonic properties of graphene and graphite from a generalized third-order ab initio approach based on the single mode relaxation time approximation and found the calculated out-of-plane k of room-temperature infinite graphite agrees with experiments, however, the in-plane k (~ 500 W/m-K) is much lower than measured, showing the single mode relaxation time approximation does not present an appropriate description of phonon thermal transport in MLG. Using an iterative solution of the linearized

ACS Paragon Plus Environment

3

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 19

BPE with harmonic and anharmonic IFCs from density functional perturbation theory (DFPT), Fugallo et al.6 found the room-temperature in-plane k converges to 2500 W/m-K at sample size of about 1000 µm for 2-layer graphene and 2200 W/m-K at sample size of about 100 µm for graphite. Their calculations involving the boundary scattering overestimated the extent of the diffusive transport regime governed by intrinsic scattering by an order of magnitude at least, as shown in this work. Despite the excellent work mentioned above, two fundamental problems concerning thermal transport in MLG remain unclear: (1) the transitions of the intrinsic in-plane k corresponding to infinite size (no scattering from boundaries) and the extent of the diffusive transport regime as n goes from n = 1 to the graphite limit and (2) tensile effects on the intrinsic in-plane k of MLG. It was reported that k of room-temperature infinite mono-layer graphene diverges under tensile strain by first principles calculations18 and MD simulations15 using the optimized Tersoff potential, however, effects of strain on MLG are still unknown. In this work, we use the selfconsistent charge density functional tight binding (SCC-DFTB) method19 to calculate the harmonic and anharmonic IFCs that govern thermal transport. We show that this method gives good accuracy and is much less computationally costly than DFPT and DFT based approaches. Then, based on iterative solutions of the linearized BPE for phonon transport, we characterize the room-temperature lattice thermal conductivities k of graphene from n = 1 to n = 4 and graphite with AB stacking at different isotropic tensile strains. We find that weak interlayer coupling decreases the extent of the diffusive transport regime by nearly two orders of magnitude and the intrinsic in-plane k by half as n goes from n = 1 to n = 2. Interestingly, there is a peak enhancement with increasing strain of the intrinsic k for the AB-stacked MLG considered and

ACS Paragon Plus Environment

4

Page 5 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

graphite, the same behavior is also observed in AA-stacked multi-layer hexagonal boron nitride (MLBN). The corresponding mechanism for the peak enhancement is discussed. For k = k αα =

the

crystallographic

1 kBT 2VN

direction

α

,

the

lattice

k

is

calculated

by

∑λ fλ ( fλ + 1)(hωλ ) vλα vλατ λ 20, 21 where ωλ , vλα and τ λ are the angular frequency, 2

the group velocity and the phonon transport lifetime, respectively, of the phonon mode λ (wavevector and branch index). The in-plane mean free path (MFP) of a phonon in mode λ is

r r defined as MFP = vλ ,in τ λ where vλ ,in is the in-plane velocity (for graphite the out-of-plane MFP r = vλ ,out τ λ ). h and fλ are the reduced Plank constant and the Bose-Einstein phonon distribution, respectively. V is the volume of the unit cell and T is temperature. For graphite, the full 3D first Brillouin zone is discretized into a Γ-centered regular grid N = N1 × N1 × N2 with N1 for the inplane directions and N2 for the out-of-plane direction (N2 > 1). For MLG and MLBN, the first Brillouin zone is 2D (N2 = 1). ωλ and vλα can be obtained from the phonon dispersions built from the harmonic IFCs, while τ λ are extracted from an iterative solution of the BPE built from the harmonic and anharmonic IFCs20,

21

considering only the intrinsic three-phonon scattering.

Technical details for the calculations of dispersion relations and three-phonon scattering rates can be found in20, 21. The harmonic and anharmonic IFCs are calculated by the SCC-DFTB method using the DFTB+ package19 with the potential derived by Gaus et al.22 and a selfconsistent charge tolerance of 1.0×10-8. A dispersion correction of UFF type as implemented in the DFTB+ package is employed to account for van der Waals interlayer interactions. An 11×11 supercell for the IFC calculations is used for mono-layer, 2-layer and 3-layer graphene. A 9×9 supercell is used for 4-layer graphene. For graphite, the corresponding supercells are 9×9×2 and

ACS Paragon Plus Environment

5

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 19

6×6×2 for the harmonic and anharmonic IFC calculations, respectively. For all cases a 1-nm cutoff is used to truncate the nearest neighbor interactions in the anharmonic IFC calculations while no additional cutoff is used in the harmonic IFC calculations. Two 1.5-nm-thick vacuum layers are used to sandwich the graphene systems to avoid periodicity effects. Before calculating IFCs, each system considered in this work is relaxed with sufficiently high precision. As an example, we obtain the equilibrium lattice constants 0.6694 nm (out-of-plane) and 0.249 nm (inplane) for unstrained graphite using a k-mesh of 21×21×5, a self-consistent charge tolerance of 1.0×10-8 and a force convergence precision of 1.0×10-6 eV/Angstrom in the DFTB+ package. These lattice parameters are in excellent agreement with corresponding experimental values 0.6711 nm and 0.2464 nm23. To be consistent with graphite, the thickness of MLG is taken as n×0.335 nm. Isotropic tensile strain, ε = ( a − a0 ) / a0 , is applied by increasing the in-plane lattice constant, a, from its equilibrium value, a0. Technically, the interlayer distance may change with strain as well, so we calculated dynamical and structural properties of graphite with ε = 0.1 for two cases: i) fixing the interlayer distance at 0.6694 nm and ii) fully relaxing the structure along the out-of-plane direction. The fully relaxed structure gives a very similar interlayer distance of 0.6675 nm and nearly identical phonon dispersions compared to the fixed case. Thus, we keep the interlayer distances fixed at 0.6694 nm for all of the graphene-based systems considered here. The calculated dispersion relations of graphite agree well with previous measurements,23, 24 as shown in Figure 1. The frequency of the high-lying flexural mode at the Γ-point is 3.55 THz, falling within experimentally observed values: 3.02 THz23 and 3.80 THz.25 As a further check of the potential and this method we compare the calculated dispersion relations for diamond with measured values,26 and obtain good agreement (see supplementary Figure S1). For unstrained graphite, we obtain approximately linear low-lying ZA dispersion in the long-wave limit,

ACS Paragon Plus Environment

6

Page 7 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

consistent with previous findings.27 For mono-layer graphene, our calculations reproduce well the quadratic ZA dispersion in the long-wave limit, while its ZA dispersion is linearized by strain, consistent with the prediction by continuum thin plate vibration theory.28 Considering the non-negligible interlayer shear coupling in AB-stacked MLG,29 the low-lying ZA modes are not decoupled from the in-plane modes. Therefore, the low-lying ZA dispersion of unstrained MLG in the long-wave limit is no longer purely quadratic. Loss of rigid rotation symmetry in ABstacked MLG also justifies this fact.

Figure 1. Comparison of calculated and measured in-plane dispersion relations for unstrained graphite. ZA1 and ZA2 denote the low-lying and high-lying flexural modes, respectively. TA and LA denote the in-plane transversal and longitudinal acoustic modes, respectively. The other branches are optical modes. For each infinite unstrained and strained system considered, the convergence of k is tested as a function of N1 (and N2 for graphite). For each N1 (N2), an iteration precision of 1×10-5 (difference of k values for successive iterative steps) is taken to ensure full self-consistent convergence of k (see supplementary Figure S2 for iteration processes in unstrained and strained 2-layer graphene with ε = 0.01 at N1 = 121).

ACS Paragon Plus Environment

7

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 19

Here we discuss thermal transport in the unstrained systems. The layer number (n) dependent intrinsic in-plane k from our calculations is shown in Figure 2a. Representatively, the convergence of unstrained 2-layer graphene as a function of N1 is shown in the inset of Figure 2a. As seen in Figure 2a, k drops by 50%, from 5500 W/m-K to 2748 W/m-K, as n goes from n = 1 to n = 2 but gradually decreases to 2350 W/m-K for n = 4, this value is very close to the calculated in-plane k = 2294 W/m-K for graphite. As mentioned by Lindsay et al.4 and Singh et al.,7 this decreasing behavior of k is due to weak interlayer coupling that breaks a phononphonon scattering selection rule in graphene and opens more scattering channels for acoustic phonons, flexural phonons in particular. Including

natural isotope (C13) scattering via

perturbation theory for a random isotope distribution,20 we obtain in-plane k = 2070 W/m-K and out-of-plane k = 7.3 W/m-K for graphite, in good agreement with experiments.10 At different temperatures (ranging from 300 K to 2000 K), we also obtain good agreement of the calculated k with measured values10 (see supplementary Figure S3). For mono-layer graphene, the gridconverged k = 5500 W/m-K is higher than recently reported k = 4300 W/m-K corresponding to 1-mm sample size by Fugallo et al..6 The boundary scattering involved in their calculations may give lower k and 1-mm sample size is not enough to give full convergence to the intrinsic k, as demonstrated by our calculations on the extent of the diffusive transport regime below. To examine this more fully, we calculate the accumulated k given as the sum of contributions from phonons with MFP less than a particular length.21 We define the saturated MFP as the length for which all mode contributions have been considered and the accumulated k is equal to the intrinsic k. Figure 2b shows the saturated MFP, representing the extent of the regime for which phonon transport is purely diffusive corresponding to the intrinsic k of the systems considered.

ACS Paragon Plus Environment

8

Page 9 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

The extent of the diffusive transport regime drops nearly two orders of magnitude, from a few centimeters to about 50 µm as n goes from n = 1 to n = 2 but gradually decreases to about 10 µm for n = 4. This is slightly larger than the value of about 8 µm for graphite. The drastic decrease in both the intrinsic k and the extent of the diffusive transport regime as n goes from n = 1 to n = 2 and their continuous transition as n goes from n = 2 to the graphite limit are consistent with the drastic dimensional transition from 2D to quasi-2D and the continuous transition from quasi-2D to 3D by increasing the layer number to infinity. The decrease of MFP results from the decrease of MFP of flexural acoustic modes while MFP (≤ 4 µm) of in-plane acoustic modes is nearly unchanged. Since the MLG sample lengths used in previous experiments are less than 6 µm,1, 2, 5, 9

the significant size effect of k in measured data may be reasonably expected. Anomalous size

dependence of k of graphene ribbons was also observed in the work by Nika et al..30 Such a size effect may explain why measured k are much lower than that of graphite for n = 4 and n ≈ 8 in experiments by Ghosh et al.1 where sample lengths of 1 ~ 5 µm were used, smaller than the extent of diffusive transport regime of graphite (~ 8 µm). We also show the temperature effect on the saturated in-plane MFP of graphite in the inset of Figure 2b. As expected, increasing temperature decreases the saturated in-plane MFP due to enhanced phonon-phonon scattering. Correspondingly, the saturated out-of-plane MFP also decreases from about 2 µm (300 K) to a few hundreds of nanometers. Further analysis shows that, although the interlayer coupling decreases k with the increase in n, the contribution to k from all the flexural phonons, denoted as n

∑k

ZAi

, is dominant (above 80%) for all the systems considered due to the large density of states

i =1

of these phonons. Moreover, the ZA branch contributions (kZAi) decrease monotonically with

ACS Paragon Plus Environment

9

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 19

increasing Γ-point frequency, i.e., kZAi > kZAj for j > i. These findings of flexural mode contributions are qualitatively consistent with those by Lindsay et al.4 and Singh et al..7 The finite sizes used in their calculations may have hindered their findings on transitions of the intrinsic k and the extent of the diffusive transport regime as n goes from 1 to infinity, as we present here.

Figure 2. (a) Dimensional transition of intrinsic k at room temperature. The inset shows k convergence for 2-layer unstrained graphene as a function of grid density N1. (b) Dimensional transition of the extent of the diffusive transport regime, i.e., the saturated in-plane MFP at room temperature. The inset shows the saturated in-plane MFP of graphite for different temperatures. We now compare the calculated k of finite 3-layer and mono-layer graphene with experimental data from Li’s group2, 8 using direct thermal bridge measurements. To be consistent with experiments, scattering mechanisms including boundary scattering, natural isotope scattering and phonon-phonon scattering are considered. The direction, x , of the temperature gradient is assumed to be along the sample length, L , and the corresponding in-plane thermal conductivity along x is defined as k = k αα . Anisotropy due to finite sample size is ignored

ACS Paragon Plus Environment

10

Page 11 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

f

because of the relatively large sizes considered21. The lifetime τ λ of a phonon in a finite sample 1

may be calculated using the Matthiessen rule4, which is expressed as

Here

1

τλ

is the intrinsic phonon-phonon scattering rate;

1

τ λiso

=

f

τλ

1

τλ

+

1 iso

τλ

+

1 b

τλ

+

1

τ λw

.

represents the scattering rate from 1

naturally occurring carbon isotopes (98.9% C12 and 1.1% C13)20, 31;

rate by contact boundaries and is expressed empirically as

respectively for a rectangular sample and a circular sample; and

1

τ λw

represents the scattering

τ λb

1

τ λb

2 vλx

=

21

L

and

1

τ λb

=

vλx

6

D

represents the scattering rate

32

due to the finite width W of the rectangular sample and is expressed as

1

τ λw

=

2 v λy W

, where v λy

is the group velocity along the width direction. We assume scattering from sample boundaries is purely diffusive. For two 3-layer graphene samples2 with widths W = 5.04 µm and W = 12.8 µm and the same length L = 5 µm, our calculations give k1 = 1853.1 W/m-K and k2 = 1951.3 W/mK, respectively. The corresponding measured values are k1 =1400±140 W/m-K and k2 = 1490±150 W/m-K. The overestimation of the calculated k values is not surprising as the added extrinsic thermal boundary resistance from the heat sensor is not considered in our calculations. Comparison between calculated and measured8 values for different lengths of mono-layer graphene at room temperature shows excellent agreement with measured data, as seen in Figure 3a. As shown in Figure 3b, we also find excellent agreement over a wide temperature range with the measured data of Chen et al.33 for mono-layer graphene suspended over holes with D=2.8 µm with natural isotopic composition. The calculated enhancement of k for isotopically purified mono-layer graphene (~12% at room temperature) is significantly below the measured values,

ACS Paragon Plus Environment

11

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 19

though within the experimental uncertainties. We note that previous calculations6,21 of k for isotopically purified mono-layer graphene also found weak enhancements similar to our results .

Figure 3. (a) Comparison of calculated k with the experimental data for different sample lengths of mono-layer graphene at room temperature. These samples have the same width W = 1.5 µm. (b) Comparison of calculated k versus temperature with measured data for mono-layer graphene suspended over holes of fixed diameter D = 2.8 µm with different isotope abundances. To investigate strain effects on the thermal transport of the systems considered here, a series of isotropic tensile strains ε are applied. For strained mono-layer graphene, we find the calculated k increases almost linearly with increasing N1 as shown in Figure 4a, thus demonstrating the intrinsic k diverges with increasing system size. This compares well with the results of k divergence of strained mono-layer graphene in previous work,15, 18 however contradicts the work of Fugallo et al.6 for which calculations using a full BPE solution gave converged intrinsic k for strained systems. Fugallo et al. used a fixed Gaussian smearing width and N1=128 in their integrations, while we used a self-adaptive Gaussian smearing method20 and N1=501 and did not find convergence with increasing N1. These differences may be the cause of the different strained k behaviors in the separate calculations for mono-layer graphene. Isotropic in-plane tensile strain

ACS Paragon Plus Environment

12

Page 13 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

does not break reflection symmetry of the mono-layer graphene system, thus the symmetry based selection rule34 for ZA modes is still applicable though the ZA dispersion is linearized.18 This reflection symmetry may be broken by irregular lattice deformation, presence of a polymer residue, or a supporting substrate. Consequently, breaking of the selection rule will open more scattering channels and reduce k, especially the ZA phonon contributions. In contrast to the mono-layer graphene, the interlayer coupling in strained MLG and graphite without the reflection symmetry opens more scattering channels and their intrinsic k remain convergent, as demonstrated in Figure 4a for 2-layer graphene, the asymptotic convergence can be seen. We note that Gu and Yang35 used the k~1/N1 convergence plot to argue that k diverges for 2D silicene. Following this approach, we achieve grid converged k of strained multi-layer graphene and graphite for sufficiently high grid densities. Similar to those in unstrained structures, flexural acoustic phonons provide the dominant contributions to k. Figure 4b shows the reduced thermal conductivity kstrained/kunstrained for 2-layer graphene, 3-layer graphene and graphite under different strains. There is a peak enhancement of ~40% corresponding to ε = 0.03 for 2-layer graphene and graphite and a peak enhancement of ~30% corresponding to ε = 0.05 for 3-layer graphene. The tensile-strain-induced improvement of intrinsic k is in contrast to those reported for other carbon-based materials such as 2D graphene nanoribbons36 and 1D carbon nanotubes,37,

38

wherein tensile strains reduce k through phonon softening mechanisms.38 For 3D diamond, our calculations give room-temperature intrinsic k = 3820.1 W/m-K and applying compressive strain enhances its k, consistent with previous DFPT(DFT) calculations.39 Moreover, applying tensile strain in diamond decreases its k significantly (see supplementary Table S1) due to softening mechanisms (see supplementary Figure S4). To ascertain the underlying mechanisms for the peak enhancements in MLG and graphite we plot lifetimes of ZA phonons with respect to

ACS Paragon Plus Environment

13

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 19

frequency and the total heat capacities of ZA modes under different strains in Figure 4c and Figure 4d, respectively, for strained 2-layer graphene. Significant enhancement in the lifetimes

Figure 4. (a) Convergence of strained mono-layer graphene and 2-layer graphene. (b) Strain dependent reduced k , i.e., the ratio of the intrinsic k under given strain to the intrinsic k without strain for 2-layer graphene, 3-layer graphene, graphite and 2-layer h-BN. Lines are presented to guide the eye. (c) Strain dependent lifetimes of low-lying ZA phonon modes for 2-layer graphene. (d) Strain dependent heat capacity of low-lying ZA modes for 2-layer graphene. of ZA phonons can be seen for ε ≤ 0.03. Such enhancement becomes smaller with increasing ε for ε > 0.03 while the total heat capacity still decreases significantly in this range, demonstrating the peak at ε = 0.03 results from the competition between the increasing lifetimes and the decreasing mode heat capacities of ZA phonons. Two points contribute to the lifetime enhancement: i) the decrease of the magnitude of the anharmonic IFCs with tensile strain21

ACS Paragon Plus Environment

14

Page 15 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

reduces the scattering matrix elements; ii) The reduction in ZA phonon density of states caused by zone center dispersion linearization18 leads to less scatterings of ZA phonons.15 The reduction in ZA phonon density of states also causes the decrease of the mode heat capacity. For other high-lying flexural modes, the variations of lifetime and mode heat capacity are similar to those of ZA modes. Our data analysis shows that the strain dependence of the group velocities of ZA phonons plays a minor role in determining the present k behavior. After exploring the peak enhancement problem in MLG, a question naturally arises: does MLBN have the same peak enhancement behavior as it is also assembled of graphene-like monolayer h-BN? Mono-layer h-BN also has quadratic dispersion relations for ZA modes.40 In order to ascertain this answer, we calculate intrinsic k of mono-layer h-BN and the reduced k of 2-layer h-BN systems at room temperature under different strains using the potential by Lukose et al..41 To validate the potential, we also calculate phonon dispersion relations of unstrained mono-layer h-BN. The same supercell and cutoff setups for IFC calculations with those in the corresponding

Figure 5. (a) Comparison of calculated and measured dispersion relations for unstrained monolayer h-BN. (b) Convergence of intrinsic k with respect to grid density for unstrained and strained mono-layer and 2-layer h-BN.

ACS Paragon Plus Environment

15

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 19

graphene systems of the same layer number are used. As shown in Figure 5a, the calculated dispersions agree with experimental values.42, 43 The infinite unstrained isotopically pure monolayer h-BN and 2-layer h-BN have grid-converged room-temperature k = 1200 W/m-K and k = 580 W/m-K, respectively, with an experimentally equivalent thickness 0.331 nm44 used for a single layer, as seen in Figure 5b. The calculated k of unstrained mono-layer h-BN agrees well with recent DFPT calculations by Cepellotti et al.,44 demonstrating the quality of this potential. Similar to that in strained mono-layer graphene, the intrinsic three-phonon scattering can not confine the intrinsic k of strained mono-layer h-BN to converge, as seen in Figure 5b for the divergent k at nonzero strains. For these cases higher order scattering (e.g., from four phonon processes) may confine the intrinsic k to converge. The k convergence of strained mono-layer graphene and mono-layer h-BN is an open question and beyond the range of this work. Representatively, the k convergence of strained 2-layer h-BN at ε = 0.01 as a function of N1 is also shown in Figure 5b for clarification. For comparison, the reduced k results for 2-layer h-BN at different strains are plotted in Figure 4b. They present similar trends to those in MLG. Therefore, we expect that the peak enhancement behavior may be an inherent quality of all multi-layer systems assembled from building blocks of purely 2D atomic layers with quadratic phonon flexural dispersions. We note that, compared with the large graphene fracture strain ~25%45 and the maximum ~40% enhancement amplitude, the almost same fracture strain46 and the ~80% enhancement amplitude observed in 2-layer h-BN indicate that engineering k of MLBN by tensile strain may be promising. In summary, from lattice dynamics calculations based on the SCC-DFTB method, we present transitions of the intrinsic k and the extent of the diffusive transport regime from 2D graphene to 3D graphite at room temperature. This knowledge enhances our understanding of dimensional

ACS Paragon Plus Environment

16

Page 17 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

and size effects on thermal transport of MLG and may guide experimental measurements. Good agreement with experimental k for graphite and finite mono-layer and 3-layer graphene samples are achieved, validating the present SCC-DFTB based calculation approach for C and BN systems and providing motivation to apply this efficient method to thermal transport calculations in other systems for which quality, well-tested potential parameters are available. Increasing tensile strain results in an interesting intrinsic in-plane k behavior in MLG and graphite due to the competition between the decreased mode capacity and increased lifetimes of flexural phonon modes, paving the way to engineer their k for thermal management applications. We also observe similar k behavior in 2-layer h-BN, indicating it may be an inherent thermal transport property for multi-layer systems assembled of purely 2D atomic layers. Further work with other MLBN and similar layered structures is expected to present a more general conclusion. ASSOCIATED CONTENT

Supporting Information This file presents i) the accuracy of the calculated phonon dispersion relations of diamond and technical details for harmonic and anharmonic IFC calculations, ii) the self-consistent iterative process of the thermal conductivity of unstrained and strained 2-layer graphene from single mode relaxation time approximation solution to full solution, iii) the accuracy of the calculated thermal conductivity of graphite under different temperatures, iv) the phonon softening as shown by phonon dispersion variations in diamond due to tensile strains and v) the calculated thermal conductivities of diamond under compressive and tensile strains. This material is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION

Corresponding Author

ACS Paragon Plus Environment

17

Nano Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 19

* E-mail: [email protected]. ACKNOWLEDGMENT We are thankful for the financial support from the Hong Kong General Research Fund under Grant Nos. 623212 and 613413. L. L. acknowledges support from the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division for work done at ORNL.

REFERENCES (1) Ghosh, S.; Bao, W. Z.; Nika, D. L.; Subrina, S.; Pokatilov, E. P.; Lau, C. N.; Balandin, A. A. Nat. Mater. 2010, 9, 555-558. (2) Wang, J. Y.; Zhu, L. Y.; Chen, J.; Li, B. W.; Thong, J. T. L. Adv. Mater. 2013, 25, 6884-6888. (3) Shahil, K. M.; Balandin, A. A. Nano Lett. 2012, 12, 861-687. (4) Lindsay, L.; Broido, D.; Mingo, N. Phys. Rev. B 2011, 83, 235428. (5) Wang, Z. Q.; Xie, R. G.; Bui, C. T.; Liu, D.; Ni, X. X.; Li, B. W.; Thong, J. T. L. Nano Lett. 2011, 11, 113-118. (6) Fugallo, G.; Cepellotti, A.; Paulatto, L.; Lazzeri, M.; Marzari, N.; Mauri, F. Nano lett. 2014, 14, 6109-6114. (7) Singh, D.; Murthy, J. Y.; Fisher, T. S. J. Appl. Phys. 2011, 110, 3622300. (8) Xu, X. F.; Pereira, L. F. C.; Wang, Y.; Wu, J.; Zhang, K. W.; Zhao, X. M.; Bae, S.; Bui, C. T.; Xie, R. G.; Thong, J. T. L.; Hong, B. H.; Loh, K. P.; Donadio, D.; Li, B. W.; Ozyilmaz, B. Nat. Commun. 2014, 5, Ncomms4689. (9) Pettes, M. T.; Jo, I. S.; Yao, Z.; Shi, L. Nano Lett. 2011, 11, 1195-1200. (10) C. Y. Ho, R. W. P., P. E. Liley., Thermal Conductivity of the Elements: A Comprehensive Review. American Chemical Society: Washington,DC, 1974. (11) Li, H.; Ying, H.; Chen, X.; Nika, D. L.; Cocemasov, A. I.; Cai, W.; Balandin, A. A.; Chen, S. Nanoscale 2014, 6, 13402-13408. (12) Zhong, W. R.; Zhang, M. P.; Ai, B. Q.; Zheng, D. Q. Appl. Phys. Lett. 2011, 98, 113107. (13) Brenner, D. W. Phys. Rev. B 1990, 42, 9458-9471. (14) Lindsay, L.; Broido, D. Phys. Rev. B . 2010, 81, 205441. (15) Pereira, L. F. C.; Donadio, D. Phys. Rev. B 2013, 87, 125424. (16) Balandin, A. A.; Ghosh, S.; Bao, W. Z.; Calizo, I.; Teweldebrhan, D.; Miao, F.; Lau, C. N. Nano Lett. 2008, 8, 902-907. (17) Paulatto, L.; Mauri, F.; Lazzeri, M. Phys. Rev. B 2013, 87, 214303. (18) Bonini, N.; Garg, J.; Marzari, N. Nano Lett. 2012, 12, 2673-2678. (19) Dolgonos, G.; Aradi, B.; Moreira, N. H.; Frauenheim, T. J. Chem. Theory Comput. 2009, 6, 266278. (20) Li, W.; Carrete, J.; A. Katcho, N.; Mingo, N. Comput. Phys. Commun. 2014, 185, 1747-1758. (21) Lindsay, L.; Li, W.; Carrete, J.; Mingo, N.; Broido, D. A.; Reinecke, T. L. Phys. Rev. B 2014, 89, 155426. (22) Gaus, M.; Goez, A.; Elstner, M. J. Chem. Theory Comput. 2012, 9, 338-354. (23) Mohr, M.; Maultzsch, J.; Dobardžić, E.; Reich, S.; Milošević, I.; Damnjanović, M.; Bosak, A.; Krisch, M.; Thomsen, C. Phys. Rev. B 2007, 76, 035439.

ACS Paragon Plus Environment

18

Page 19 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nano Letters

(24) Maultzsch, J.; Reich, S.; Thomsen, C.; Requardt, H.; Ordejón, P. Phys. Rev. Lett. 2004, 92, 075501. (25) Nicklow, R.; Wakabayashi, N.; Smith, H. Phys. Rev. B 1972, 5, 4951-4962. (26) Schwoerer-Böhning, M.; Macrander, A.; Arms, D. Phys. Rev. Lett. 1998, 80, 5572-5575. (27) Neto, A. C.; Guinea, F.; Peres, N.; Novoselov, K. S.; Geim, A. K. Rev. Mod. Phys. 2009, 81, 109-162. (28) Virgin, L. N., Vibration of axially-loaded structures. Cambridge University Press: 2007. (29) He, X.; Kitipornchai, S.; Liew, K. Nanotechnology 2005, 16, 2086-2091. (30) Nika, D. L.; Askerov, A. S.; Balandin, A. A. Nano Lett. 2012, 12, 3238-3244. (31) Lindsay, L.; Broido, D. A.; Mingo, N. Phys. Rev. B 2011, 83, 235428. (32) Turney, J.; McGaughey, A.; Amon, C. J. Appl. Phys. 2010, 107, 024317. (33) Chen, S. S.; Wu, Q. Z.; Mishra, C.; Kang, J. Y.; Zhang, H. J.; Cho, K. J.; Cai, W. W.; Balandin, A. A.; Ruoff, R. S. Nat. Mater. 2012, 11, 203-207. (34) Lindsay, L.; Broido, D. A.; Mingo, N. Phys. Rev. B 2010, 82, 161402. (35) Gu, X.; Yang, R. J. Appl. Phys. 2015, 117, 025102. (36) Wei, N.; Xu, L.; Wang, H.-Q.; Zheng, J.-C. Nanotechnology 2011, 22, 105705. (37) Li, X. B.; Maute, K.; Dunn, M. L.; Yang, R. G. Phys. Rev. B 2010, 81, 245318. (38) Xu, Z. P.; Buehler, M. J. Nanotechnology 2009, 20, 185701. (39) Broido, D. A.; Lindsay, L.; Ward, A. Phys. Rev. B 2012, 86, 115203. (40) Kern, G.; Kresse, G.; Hafner, J. Phys. Rev. B 1999, 59, 8551-8559. (41) Lukose, B.; Kuc, A.; Frenzel, J.; Heine, T. Beilstein J. Nanotechnol. 2010, 1, 60-70. (42) Nemanich, R.; Solin, S.; Martin, R. M. Phys. Rev. B 1981, 23, 6348-6356. (43) Rokuta, E.; Hasegawa, Y.; Suzuki, K.; Gamou, Y.; Oshima, C.; Nagashima, A. Phys. Rev. Lett. 1997, 79, 4609-4612. (44) Cepellotti, A.; Fugallo, G.; Paulatto, L.; Lazzeri, M.; Mauri, F.; Marzari, N. Nat. Commun. 2015, 6, Ncomms7400. (45) Zhang, J.; Zhao, J.; Lu, J. Acs Nano 2012, 6, 2704-2711. (46) Peng, Q.; Ji, W.; De, S. Comput. Mater. Sci. 2012, 56, 11-17.

Table of Contents graphic

ACS Paragon Plus Environment

19