Parametrization of Density Functional Tight-Binding Method for

Subscriber access provided by Warwick University Library

Article

Parameterization of Density Functional Tight Binding Method for Thermal Transport in Bulk and Low-Dimensional Si Systems Qi Wang, Xinjiang Wang, Ruiqiang Guo, and Baoling Huang J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b04182 • Publication Date (Web): 22 Jun 2017 Downloaded from http://pubs.acs.org on June 25, 2017

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.

The Journal of Physical Chemistry C 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 49

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

The Journal of Physical Chemistry

Parameterization of Density Functional Tight Binding Method for Thermal Transport in Bulk and Low-dimensional Si Systems Qi Wang1, Xinjiang Wang1, Ruiqiang Guo1, Baoling Huang1,2*

1

Department of Mechanical and Aerospace Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong

2

The Hong Kong University of Science and Technology Shenzhen Research Institute, Shenzhen, 518057, China

ABSTRACT

Understanding the thermal transport in different Si-based materials is of both practical and academic importance due to their essential role in modern electronic, microelectromechanical applications and other areas. Conventional atomic modeling approaches such as the Density Functional Theory offer high accuracy but can hardly handle a large system while empirical potentials adopted by classical molecular dynamics often lack accuracy or transferability. We 1

ACS Paragon Plus Environment


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 49

developed a new parametrization of Si-Si interaction of the Density Functional based tight binding method for the atomic-scale investigation of thermal transport properties in various Si systems. It shows that this parameterization can accurately predict many harmonic and anharmonic thermal transport properties in different silicon systems such as single crystalline silicon, silicene and silicene nanoribbons, showing excellent computation efficiency and transferability. Therefore, the Si-Si parameter set can contribute to the fundamental understanding of thermal transport in various complex Si-based systems which are difficult to accurately model using traditional methods.

*E-mail: [email protected]

1. INTRODUCTION Silicon, the second most abundant element on earth, has been an essential component in modern transistors, sensors and actuators, and numerous new applications with bulk silicon and low-dimensional Si-based structures such as silicon nanowire, superlattice, and 2D silicene continue emerging. Understanding the thermal transport in these Si-based devices and systems is critical for their design and thermal management. Meanwhile, various Si-based systems also provide an excellent platform for investigating fundamental thermal transport mechanisms at different scales. To obtain detailed mode-wise thermal transport information, atomic-scale

2


Page 3 of 49

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


approaches such as molecular dynamics (MD) simulations and Boltzmann Transport Equation (BTE) based lattice dynamics methods are widely adopted. The accuracy of these atomic-scale methods essentially relies on the quality of the interatomic potentials, which can be either empirically fitted from experimental results1-3 or directly calculated by density functional theory (DFT)4,5. The empirical interatomic potentials are widely used by classical MD simulations and lattice dynamics calculations because they are very computationally efficient and easy to implement. However, these empirical potentials are often fitted with selected properties of a given structure without explicitly considering the atomic bonding variation among different structures, and therefore their accuracy and transferability are often serious concerns6,7. For example, even for single crystalline silicon, many classical potentials that can well predict mechanical properties of Si, e.g., Tersoff8 and Stillinger-Weber2 (SW), often give poor predictions on the phonon dispersion relations as well as the lattice thermal conductivity at the room temperature9,10. Meanwhile, most of these empirical potentials developed for single crystalline silicon cannot accurately predict the structures or the properties of Si systems involving many surface atoms such as silicene11, the silicon counterpart of graphene. In contrast, the DFT method directly considers the electron density distribution and requires no adjustable parameters, providing high accuracy and excellent transferability. It has been successfully implemented with BTE to predict the thermal transport 3



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 49

properties of many bulk materials such as diamond12 and silicon13 as well as 2D materials like planar graphene14 and buckling MoS215. However, the DFT method is very computationally demanding and currently the applications of DFT method mainly concentrate on small systems like simple crystals. The low computational efficiency of the DFT method therefore limits its applications for thermal transport property investigations of large systems and complex structures such as alloys and nanostructures. Therefore a method that provides a good balance between accuracy, speed and transferability is highly desirable for the investigations of such systems. The density functional tight binding (DFTB) method, as an approximation to the DFT method, has been proved successful in chemistry, physics and biology. It treats the electron density as a function of predetermined atomic orbitals, whose parameters and couplings are pre-calculated and tabulated. Therefore, it requires no integral evaluation during runtime and is almost two orders of magnitude faster than the conventional DFT method. It can be an efficient quantum simulation choice for the thermal property calculations of large systems or complicated structures with a high computational efficiency. Due to its origin from DFT, the reasonable accuracy and transferability could also be expected with suitable Slater-Koster files (namely parameters). Several packages incorporating the DFTB method such as DFTB+16 and hotbit17 have been served for years for various purposes along with many sets of parameters. More details about 4


Page 5 of 49

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


DFTB applications can be found elsewhere17,18 . However, up until now, only a few lattice thermal conductivity calculations have been conducted with the DFTB method19,20. What makes it difficult is that the DFTB parameterization is often conducted with limited fitting data sets and in special purposes such as achieving accurate electronic parameters or surface energy calculations, and thermal transport properties are rarely considered in the past parameterizations. For example, as will be shown later, the Si-Si interaction parameters provided by DFTB+ for solids and surfaces (denoted as “pbc” set21) can accurately reproduce the lattice structure and mechanical properties of single crystalline silicon but significantly overestimate its thermal conductivity at 300 K (302 W/m-K) compared with the experimental value22 (150 W/m-K). Therefore such parameter sets are unsuitable for the quantitative study of thermal transport properties in various Si systems, and new parametrization for Si-Si interaction with good accuracy and transferability is needed. In this study, we propose a new parametrization for the Si-Si interaction within the DFTB scheme for the study of thermal transport properties in different Si systems. We have extensively tested this parametrization of Si-Si interaction (named as “SiTherm” parameters) for silicon, 2D silicene and silicon nanoribbons (SiNRs), and the SiTherm parameters are proved to be able to accurately predict the lattice structures and many thermal transport properties of various Si systems. This new set of Si-Si parameters demonstrates excellent accuracy and transferability, and 5



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 49

the computation speed of DFTB with the SiTherm set is also much faster than the conventional DFT method (about 2500 times faster for silicene systems and 20 times faster for silicon systems), rendering it an ideal option for the study of large-scale Si systems. 2. THEORY

2.1 Density Functional Based Tight Binding Method The self consistent charge density functional based tight binding (SCC-DFTB) total energy can be regarded as the second-order Taylor expansion of the DFT total energy around a reference

r density n0 (r ) . This method treats the electron density as a function of predetermined atomic orbitals. The total energy expression in the DFTB scheme can then be expressed as23

Etotal = ∑ i

occ

Ψi H0 Ψi +

1 N ∑ γ αβ ∆qα ∆qβ + Erep 2 α ,β

(1)

where Ψ i is the Kohn-Sham orbital, H 0 is the normal Hamilton operator depending only on

r n0 (r ) , γ αβ represents a function of distance Rαβ and Hubbard parameter U, and ∆q is the Mulliken charge. The remaining term

Erep

stands for the ion-ion interaction and

exchange-correlation interactions, which can be roughly referred as the summation of repulsive r interaction Eαβ between the atoms α and β , i.e.,

6


Page 7 of 49

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


Erep =

1 ∑ Eαβr ( Rαβ ) 2 αβ

(2)

Introducing the minimal local basis of atomic orbitals η µ to expand the Kohn-Sham orbitals ( Ψi = ∑ cµi ηµ , µ = s, p, d and f), and combining the expansion of repulsive part, the SCC-DFTB total energy becomes:

Etotal = ∑i

occ

1

1

cµ cν H µν + ∑α β γ αβ ∆qα qβ + ∑αβ Eαβ ( Rαβ ) . ∑∑ 2 2 µ ν i

i

0

N

r

(3)

,

0 Here H µν is related to the Kohn-Sham equations in DFTB:

αβ 0 H µν = H µν +

1 ∑ Sµν ∆qk (γ α k + γ β k ), µ ∈ α ν ∈ β . 2 k

(4)

When introducing the localized atomic orbitals, due to the variational principle their coefficients cνi should be determined by:

∑ν cνα ( H αβ − ε α S µν ) = 0 ,

(5)

where S µν represents the overlap matrix element. Equations (4) and (5) in general are solved iteratively to ensure self-consistency. 2.2 SCC-DFTB parametrization

7



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 49

According to the SCC-DFTB theory, there are four traditional key steps when parametrizing a pair of interactions24. These steps involve 1) the DFT calculation of the localized basis functions

r nµ and the reference density n0 (r ) , 2) the determination of confining potential Vconf (which is used for the determination of localized atomic orbitals), 3) tabulating the Hamiltonian H µν and overlap matrix elements S µν and 4) fitting the repulsive potential Erep as the spline functions or polynomials. From a more practical standpoint, the confining potential Vconf and the repulsive potential Erep are important and easy to adjust referring to the starting goal, for example, phonon

dispersion relation. Vconf

r usually takes the form of    r0 

N

where N = 2 is mostly approximated

and r0 is chosen to be the double of covalent radius empirically. As for the term Erep , in the fitting process it always performs as the difference of DFT total energy and the SCC-DFTB electronic energy (the first two terms in Equation (1) on the right side) for a reference structure:

ele Erep = EDFT (r ) − ESCC − DFTB (r ) .

(6)

r ( Rαβ ) } could be collected by Thus, for a given equilibrium structure, data points {R, Eαβ

stretching the bond between the atoms α and β near the equilibrium length for the energy curve fitting with polynomials. More details about DFTB and the parametrization can be found elsewhere17,18,24-26. 8


Page 9 of 49

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


3. COMPUTATIONAL DETAILS

To test the accuracy and transferability of the newly developed DFTB parameter set SiTherm, we have calculated various structural, harmonic and anharmonic lattice properties, and thermal conductivities of different Si systems including single crystalline silicon, 2D silicene and silicene nanoribbons using the DFTB method with the SiTherm set. These results were compared with those from experiments, DFT, empirical force fields and other common DFTB parameter sets (e.g., pbc and matsci27) in the literature. The detailed computation methods and details are introduced below. 3.1 Structure A fundamental criterion to evaluate the accuracy and transferability of a force field is its capability to accurately predict the lattice parameters of different structures, which may significantly influence the calculations of many other properties. We have optimized the structures of bulk silicon, silicene and silicene nanoribbons using the DFT and DFTB methods. For silicon structure optimization, within the DFT framework, a conventional unit cell of silicon was relaxed under both the local density approximation (LDA) of Perdew and Zunger28, and the Perdew-Burke-Ernzerhof (PBE)29 form of the generalized gradient approximation (GGA) for the projector augment wave method (PAW)30 implemented in the VASP. The plane wave cutoff was 9



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 49

set to 250 eV. The equilibrium lattice parameters were obtained by fitting the energy over a reasonable large volume range according to the Birch-Murnaghan equation of state31 . Within the DFTB framework, the unit cell was relaxed using the same method with different parameter sets including the pbc, matsci and SiTherm sets. All these relaxations adopted a 25×25×25 k-point sampling grid, with a consistent electronic convergence criterion of 10-9 eV. For silicene structure optimization, within both the DFT and DFTB frameworks, a vacuum layer of a thickness of 20 Å was adopted to avoid the adjacent layers’ interactions due to the periodic boundaries of the simulation cell. Then the primitive unit cell of silicene was relaxed with a convergence criterion of 10-6 eV/Å for the forces on each atom, and the electronic convergence criterion was 10-9 eV with a 29×29×1 k-points sampling grid. In the DFT calculations the plane wave cutoff energy was set as 550 eV. Three different parameter sets including the pbc, matsci and SiTherm sets were adopted in these DFTB calculations. The SiTherm set was also applied in the structure optimizations of pristine zigzag and armchair silicene nanoribbons (ZSiNRs and ASiNRs) in the DFTB scheme. The (2×1) supercells 32 of SiNRs were considered with a force convergence criterion of 10-4 eV/Å and an electronic convergence criterion of 10-5 eV. 3.2 Cohesive energy

10


Page 11 of 49

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


Cohesive energy for a solid is the energy required to separate its constituents into free neutral atoms, usually defined as the difference of energy of a free atom and crystal energy (N is the number of atoms in a unit cell):

Ecoh = E free atom −

1 Eunit . N

(7)

The cohesive energy indicates the interatomic bonding strength and is correlated with the lattice structure that a solid adopts. We calculated the cohesive energy of single crystalline silicon using both the DFT and the DFTB methods with different parameter sets. A k-point sampling grid of 35×35×35 was chosen for the static energy calculations of one silicon conventional unit cell in our DFT and DFTB calculations. 3.3 Harmonic and anharmonic lattice properties Thermal transport is directly influenced by both the harmonic and anharmonic lattice properties, and a convenient approach to judge the quality of a force field is to see if the force field can accurately yield these lattice properties (e.g., the elastic moduli and Grüneisen parameters) due to the computational simplicity. Young’s modulus and shear modulus can be derived from the strain derivative of the total energy of a lattice33. In this study, a conventional unit cell of silicon was adopted for the calculations of lattice energy under different strains using DFTB+ and VASP. Also,

11



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 49

using the DFTB method with the SiTherm set, we calculated the force constant K and in-plane stiffness C of a (8×1) supercell of 10-ASiNR under uniaxial tension along the ribbon length direction. Here K was calculated according to K =

∂ 2 Es ( Es is the strain energy, and c is the ∂c 2

corresponding lattice constant) and the in-plane stiffness C is defined as C =

1 ∂ 2 Es (A0 is the A0 ∂ε 2

equilibrium area of the supercell, ε is the strain and given by ε = (c − c0 ) / c0 , and c0 is the lattice constant of the unstrained supercell). On the other side, the Grüneisen parameter is often regarded as the indicator of lattice anharmonicity and is correlated with the Umklapp scattering strength. The mode Grüneisen parameter is defined as γ λ = − d [ln ωλ ] / d [ ln Ω ] , where ω λ is the phonon frequency of phonon mode λ , and Ω is the unit cell volume. Therefore, the mode Grüneisen parameters could be extracted from the frequency shifts of phonon modes under varying lattice volumes (±0.5% in this study). 3.4 Lattice thermal conductivity Based on the phonon BTE and Fourier’s law, the lattice thermal conductivity tensor can be calculated by summing up the contributions from all phonon modes q, i.e.,

κ ij = −

1 ∑ hωq vqi nqo ( nqo + 1) Fqj NoΩ q

,

(8)

12


Page 13 of 49

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


where

v is

the phonon

nqo = 1 − exp ( −hω / ( k BT ) )

−1

group

velocity vqi = ∂ω ( q ) / ∂i q ,

Ω

is

unit

cell

volume,

is the equilibrium phonon Bose-Einstein occupancy, and No is the total

number of q points (within the first Brillouin zone), and F can be obtained by solving the BTE −v i ( q ) nq0 ( nq0 + 1)

hωq k BT 2

=

1 ∑ ( Pqq,q'' ' + Pqq,q' '' + Pqq',q '' )( Fqα + Fqα' + Fqα'' ) + Pqb Fqi 2 q ', q '' 

Pqq,q'' ' = 2π f qo f qo' ( f qo" + 1) V3 ( q, q ', q '' ) δ (ωq + ωq ' − ωq '' ) 2

,

Pqb =

vq L

,

(9)

nqo ( nqo' + 1)

(10)

where Pqq,q'' ' Pqq,q' '' and Pqq',q '' are three-phonon scattering probability,

Pqb

is the boundary scattering

probability, δ is the Dirac delta function, V3 is the anharmonic force constants projected onto the eigenvector space, and L is the Casimir effective length of the sample. Here we only considered the intrinsic phonon-phonon scattering and phonon-boundary scattering as an example, other scattering mechanisms can also be included in Eq.(9) through the Matthiessen rule. Equation (8) can therefore be solved iteratively using the preconditioned conjugate-gradient method34 along with an adaptive smearing scheme35. To ensure the convergence of intrinsic thermal conductivity, careful convergence tests were conducted and a 29×29×29 q mesh for silicon and a 129×129×1 q mesh for silicene were selected respectively. More technical details can be found in our previous papers19,36,37. Besides, for silicene the monolayer thickness was taken as 0.42 nm, the equilibrium van der Waals diameter of silicon atoms.

13



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 49

The only required inputs to calculate the lattice thermal conductivity are therefore the harmonic and anharmonic interatomic force constants (IFCs), which are mainly used for extracting phonon dispersion relations, group velocities and phonon scattering rates. Usually IFCs can be acquired from the first principle approaches such as DFT calculations or extracted from classical MD simulations. In our DFTB calculations, the IFCs were obtained from the SCC-DFTB calculations using DFTB+. The harmonic interatomic force constants (IFCs) are defined as:

Ψ αβ (lb, l ' b ') =

∂ 2V , ∂uα (lb )∂u β (l ' b ')

(11)

where V denotes the potential energy, uα (lb) is the α coordinate of the bth atom inside the lth unit cell. Similarly, the anharmonic IFCs are defined as:

Ψ αβγ (lb, l ' b ', l " b ") =

∂ 2V . ∂uα (lb)∂u β (l ' b ')∂uγ (l " b ")

(12)

For silicon, harmonic and anharmonic IFCs were obtained by constructing a cubic 2×2×2 supercell with a 9×9×9 k-point sampling grid, while for silicene a 9×9×1 supercell of primitive unit cell was used with a 3×3×1 k-point sampling grid to extract harmonic IFCs and a 5×5×1

14


Page 15 of 49

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


supercell with a 5×5×1 k-point sampling grid to extract anharmonic IFCs. For all the self-consistent electronic calculations, the convergence criterion was set as 10-9 eV. 4. RESULTS AND DISCUSSIONS

4.1 Silicon 4.1.1 Lattice structure, cohesive energy and elastic moduli of bulk silicon We first calculated the equilibrium lattice structure and some macroscopic mechanical properties of single crystalline silicon using the DFT and DFTB methods. Both the LDA and GGA (PBE) approximations were adopted in the DFT calculations. Two widely used DFTB parameter sets for Si, i.e., the pbc and matsci sets, were also employed here for comparison. Table 1 shows the optimized lattice constants, cohesive energy, Young’s modulus and shear modulus of single crystalline silicon predicted by the DFT and DFTB methods with different sets of parameters, along with the experimental values38-40. For lattice constants, the DFT/LDA(PBE) results are close to previous DFT calculations41, and also agree quite well with the experimental value38 5.4307 Å despite the slight underestimation from the LDA treatment and overestimation from the PBE treatment. The lattice constant calculated using the SiTherm set is in excellent agreement with the measured value38 (only 0.2% difference) while the other two DFTB parameter sets, pbc and matsci,

15



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 49

give relatively larger discrepancies compared with the measured value (0.54% and 5.2% , respectively). Cohesive energy characterizes the interatomic bonding strength and is an important parameter in the studies of surface and defect formation and propagation. We first calculated the silicon cohesive energies with DFT/PBE and DFT/LDA, which are almost the same as previous DFT results42. As shown in Table 1, the result from the DFT/PBE combination (4.60 eV) agrees the best with the experimental value (4.693 eV). The cohesive energy from the DFTB/SiTherm combination is 4.12 eV, slightly smaller than the experimental value but still better than the predictions from DFT/LDA (14% deviation) and DFTB/matsci (64% deviation). Elastic properties of solids are important in material science and engineering, and they can reveal many fundamental solid-state properties such as phonon group velocity and thermal expansion. Table 1 also compares Young’s modulus and shear modulus predicted from different methods. The calculation results from all the DFT and DFTB methods except the DFTB/matsci combination fall within the experimental data range40. DFTB predictions with the SiTherm and pbc parameter sets match quite well with the DFT results. 4.1.2 Phonon dispersion relation, DOS, and Grüneisen parameters of bulk silicon

16


Page 17 of 49

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


In lattice dynamics, the phonon dispersion relation plays an essential role in calculating the thermal conductivity. From the phonon dispersion relation, the dynamical properties of a lattice like phonon frequencies and group velocities could be extracted. Figure 1 compares the phonon dispersion relations along some high-symmetry directions obtained with the DFT method and the DFTB method (using pbc, matsci and DFTB/SiTherm sets) together with the experimental results43,44. The corresponding phonon DOS are also shown in the right panel. The accuracy of the first-principles methods have been proved in predicting the phonon dispersion relation of bulk silicon45, as also confirmed in our comparison. But for DFTB parameter sets pbc and matsci, the discrepancies between the calculated results and experimental ones appear much larger. For example, although the pbc set predicts acceptable results for optical phonons, the dispersion relation of acoustic phonons, which contributes the most to the thermal conductivity, is not well predicted. In contrast, the calculated phonon dispersion relation from the SiTherm set, especially for acoustic phonons, is much improved and is in excellent agreement with the experimental data and DFT results. Compared with the other two sets of DFTB Si-Si parameters this improvement will consequently lead to a better prediction of the silicon thermal conductivity. Such improvements can also be found in the DOS plot. The DOS predicted by DFTB/SiTherm almost overlaps with that from DFT/PBE method over most of the frequency regime. In addition, the 17



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 49

optical branches predicted by the DFTB/SiTherm combination seem to be slightly flat compared with the DFT/PBE and DFTB/pbc combinations, which may stem from the inclusion of d-orbital in the parameterization. But this is expected to have minor influence on the thermal transport properties since acoustic phonons dominate the thermal transport in bulk Si. Figure 2 shows the comparison of phonon dispersion relation curves calculated with the classical empirical potentials including the SW potential and Tersoff potentials and those from the DFTB method with the SiTherm parameter set together with the experiment data. Obviously, the widely used SW and Tersoff potentials overpredict frequencies of both acoustic modes and optical modes beyond the Γ point due to their short-range feature in the harmonic approximation, as also observed and explained by Poter et al46. Thus, given the comparable calculation speed of the DFTB method with the classical MD, property calculations related to harmonic behaviors of Si systems can recur to the DFTB method with the SiTherm set. Grüneisen parameter is a good measure of anharmonic scattering strength of phonons. Figure 3 shows the mode-wise Grüneisen parameters calculated with the DFT method and the DFTB method with the SiTherm set, and the experiment results5 are also shown. Overall, the results predicted with the DFTB/SiTherm set, especially those for acoustic phonons, agree quite well with those experimental and DFT results. Although the results for optical phonons are slightly 18


Page 19 of 49

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


overestimated (~10%), the influence is believed to be minor, especially considering the fact that optical phonons contribute little to the thermal transport in silicon. The accurate prediction of mode-wise Grüneisen parameters indicates that the SiTherm parameter set can well capture the anharmonic behaviors of phonons and is therefore suitable for the studies related to lattice anharmonicity. 4.1.3 Thermal transport in bulk silicon With all the ingredients (the optimized unit cell, the harmonic IFCs related to phonon dispersion and the anharmonic IFCs related to the scattering strength) ready, the thermal conductivities of silicon under different temperatures were then calculated. Figure 4 shows the thermal conductivities of single crystalline silicon calculated by the DFTB and DFT methods from 200 to 700 K, along with the experiment values22. As mentioned before, the DFTB/pbc combination overestimates the silicon thermal conductivity by a factor of 2 at 300 K while DFTB/matsci gives even poor results (not shown here). In contrast, DFTB/SiTherm combination can accurately predict the thermal conductivities of silicon over a wide temperature range, well matching with the DFT and experiment results. The average deviations from the experimental values and DFT results at 300 K are less than 7.1% and 2.9%, respectively. The phonon relaxation times of bulk silicon calculated by the DFT and DFTB methods are also plotted in Fig. 5. It is clear that the phonon 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

Page 20 of 49

mode relaxation times from DFTB/SiTherm are generally in good agreement with DFT/PBE results over a wide frequency range. DFTB/PBE and DFTB/SiTherm predict quite similar frequency dependences of the relaxation times for acoustic phonons ( τ

p− p

ω −3.0 and τ p − p

ω −3.2 ,

respectively) at 300 K. It can therefore be concluded that the SiTherm parameter set can be used to study the mode-wise phonon transport in Si. 4.2 2D silicene 4.2.1 Lattice structures of silicene and SiNRs Silicene is a 2D material of silicon, an analogue of graphene, except that its buckling doesn’t exist in graphene. Silicene has attracted wide interest since its successful synthesis in 201047. Recently, monolayer silicene has been applied to field-effect transistors operating at room temperature48. Its thermal conductivity has been extensively investigated using the classical MD49,50 and DFT/BTE method11,51,52. To test the transferability of the SiTherm set for the study of heat transfer in low-dimensional Si systems, relevant properties of silicene were also calculated using the DFTB method with this parameter set. Table 2 lists the lattice constants and buckling heights of the optimized silicene structure from both the DFT and DFTB methods. The results calculated by the DFT methods are similar to the literature values11,51,52. The SiTherm set can also accurately yield the silicene lattice structure, especially the buckling height, with only 2.3% and 20


Page 21 of 49

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


0.44% deviations from the DFT/LDA and DFT/PBE results. In contrast, DFTB/pbc overestimates the buckling height (more than 37% deviation from the DFT results) and DFTB/matsci virtually couldn’t predict the buckling structure. In practice, SiNRs are probably the most common forms of silicene. It is also of importance to test the transferability of the SiTherm parameters with SiNRs. As for pristine ZSiNR or ASiNR optimizations, two samples with (2×1) supercells were optimized using the DFT/LDA and DFTB/SiTherm methods, and optimized pristine SiNR structures are shown in Fig. 6. Obvious edge reconstructions can be observed in both DFT and DFTB results, in good agreement with previous DFT calculations for SiNRs32,53. However, neither the pbc nor the matsci set could reproduce the geometries of the SiNRs or the edge reconstructions. 4.2.2 Harmonic lattice properties of silicene and SiNRs For 2D materials such as graphene and silicene, the phonon dispersion relation (especially ZA branch in the long wavelength limit) plays an essential role in determining the lattice thermal conductivity54,55. To obtain accurate phonon dispersion relations and thermal conductivities of 2D materials, in our DFT and DFTB calculations, the harmonic and anharmonic IFCs were applied with symmetry constrains including the permutational, translational, and rotational invariances56. Figure 7 shows the calculated phonon dispersion relation and DOS of silicene from both DFT and 21



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 22 of 49

DFTB methods, and the inset is the logarithmic plot of the dispersion relation of low frequency phonons near the Γ point. The calculated LA and TA branches of silicene from the DFT/LDA method show a linear behavior near the Γ point, while the ZA branch shows a quadratic relation, in line with previous first-principles studies on silicene52,57 and other 2D materials 58,59. The phonon dispersions from DFTB/SiTherm are in good agreement with the DFT/LDA method over a wide frequency range, especially near the Γ point, demonstrating the ability of capturing the long wavelength phonons’ behavior. The DOS from the SiTherm set in Fig. 7 also proves the overall agreement with that of DFT/LDA. Compared with the matsci set, the phonon dispersion relation of silicene predicted by SiTherm set shows reasonable optical frequencies and no imaginary frequencies, while comparing with the pbc set, the LA branch from SiTherm set is much improved. The fundamental mechanical properties, such as the force constant and in-plane stiffness, of a pristine 10-ASiNR in the linear deformation regime were also calculated with the SiTherm set using the DFTB method. The force constant and in-plane stiffness from DFTB/SiTherm are 27.6 N/m and 53.3 N/m, similar to previous DFT values 30 N/m and 51 N/m60, which shows that the SiTherm set is applicable for Si nanosturcture mechanical property calculations in the harmonic regime. 4.2.3 Thermal transport in silicene 22


Page 23 of 49

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


The thermal conductivity of silicene was calculated using the DFT/LDA and the DFTB/SiTherm methods with the iterative method, which is necessary because of the much underestimation of 2D materials’ thermal conductivity from the single-mode relaxation time approximation method51,59. Figure 8 shows the thermal conductivities of the infinite silicene over a temperature range from 200 to 800 K from both the DFT/LDA and DFTB/SiTherm methods. Our DFT calculations with the PBE and LDA functionals produced similar results, so only the DFT/LDA results are shown in Fig. 8. Previous DFT studies11,51,52,55,57 with LDA or PBE predicted that the thermal conductivity of silicene was normally in the range of 20-36 W/mK at 300 K. The thermal conductivity of silicene predicted by DFT/LDA in our calculations is 33.3 W/mK at 300 K, agreeing well with the literature results, while that from DFTB/SiTherm method is 33.6 W/mK at 300 K, quite close to the DFT/LDA results. In contrast, the pbc set gives a value of 79 W/mK, much higher than the DFT results. The temperature-dependent thermal conductivities from both methods are also in good agreement (the largest deviation is only 8.5% at 800 K), as shown in Fig.8. Figure 9 shows the calculated phonon relaxation times of silicene at 300 K. In our DFT calculations, the relaxation times of LA and TA modes have quite weak frequency dependence while the ZA phonon modes dominate the thermal transport in silicene. It is noted that the role of the ZA modes in determining the thermal transport in silicene is still in debate

11,50,51,55

. The discrepancy may be due to the 23



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 24 of 49

different methods used (e.g., classical MD vs first-principles), different interatomic potentials, and different treatments of the dispersion near the Γ points, which will influence the behaviors of long-wavelength phonons. Recent first-principle works52,57 using the all-symmetry-enforced harmonic interatomic force constants also predicted a thermal conductivity of around 22-29 W/m-K for silicene at 300 K and they concluded that ZA modes dominate the thermal transport. Apparently, with the same treatment, our DFTB results show excellent agreements with the DFT ones. The relaxation times of low frequency ZA phonon modes follow a power law behavior ( τ ∝ ω −1.17 ), a weaker frequency dependence than that for graphene ( τ ∝ ω −1.4 )54, implying a convergent thermal conductivity for infinite sample. The higher scattering strength at the long wavelength limit is believed to be related to the buckling structure. The DFTB/SiTherm method can also accurately predict the phonon mode relaxation times over the entire the frequency range compared with the DFT/LDA method, especially in the long wavelength limit ( τ ∝ ω −1.00 ), which indicates that DFTB/SiTherm can accurately capture the mode-wise phonon transport behavior, especially for the low frequency ZA phonon modes, which play a crucial role in accurately determining lattice thermal conductivity. The accordance of mode-wise phonon relaxation times between DFTB/SiTherm method and DFT/LDA method shows the capability of the SiTherm set in capturing detailed phonon transport properties of Si-based nanostructures. Combining with the 24


Page 25 of 49

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


good performance in the geometry optimization and phonon dispersion relation calculation, the transferability of SiTherm set for predicting harmonic and anharmonic properties from the bulk to the 2D is demonstrated. 5. CONCLUSION

We have developed a new set of Si-Si DFTB parameters, SiTherm, for the study of thermal transport properties of different Si systems. Detailed tests on the structural, mechanical, harmonic and anharmonic properties of bulk silicon, 2D silicene and SiNRs have been conducted to verify the efficiency, accuracy and transferability of the proposed parameter set. Comparison with DFT, classical potentials and experimental results shows that this new parameter set can provide accurate predictions of thermal transport properties of various Si systems with a relatively low demand for computational resource. The SiTherm parameter set has been proved to be capable of describing both the Si atoms in a lattice and those on a surface, and it is expected to be applicable for other low-dimensional structures such as Si nanowires. The good balance between the accuracy, speed and transferability of this new parameter set make it an efficient choice for the study of thermal transport in large Si systems.

ACKNOWLEDGMENTS 25



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 26 of 49

We are thankful for the financial support from the National Natural Science Foundation of China (Grant No. 51376154) and the Hong Kong General Research Fund (Grant Nos. 613413 and 623212).

26


Page 27 of 49

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


REFERENCES

1. Tersoff, J. New Empirical Model for the Structural Properties of Silicon. Phys. Rev. Lett. 1986, 56, 632. 2. Stillinger, F. H.; Weber, T. A. Computer Simulation of Local Order in Condensed Phases of Silicon. Phys. Rev. B 1985, 31, 5262. 3. Daw, M. S.; Baskes, M. I. Embedded-Atom Method: Derivation and Application to Impurities, Surfaces, and Other Defects in Metals. Phys. Rev. B 1984, 29, 6443. 4. Ward, A.; Broido, D. Intrinsic Phonon Relaxation Times from First-Principles Studies of the Thermal Conductivities of Si and Ge. Phys. Rev. B 2010, 81, 085205. 5. Esfarjani, K.; Chen, G.; Stokes, H. T. Heat Transport in Silicon from First-Principles Calculations. Phys. Rev. B 2011, 84, 085204. 6. Luo, T.; Chen, G. Nanoscale Heat Transfer–from Computation to Experiment. Phys. Chem. Chem. Phys. 2013, 15, 3389-3412. 7. Huang, B.; Kaviany, M. Ab Initio and Molecular Dynamics Predictions for Electron and Phonon Transport in Bismuth Telluride. Phys. Rev. B 2008, 77, 125209. 8. Tersoff, J. New Empirical Model for the Structural Properties of Silicon. Phys. Rev. Lett. 1986, 56, 632. 9. Volz, S. G.; Chen, G. Molecular-Dynamics Simulation of Thermal Conductivity of Silicon Crystals. Phys. Rev. B 2000, 61, 2651. 10. Schelling, P. K.; Phillpot, S. R.; Keblinski, P. Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity. Phys. Rev. B 2002, 65, 144306. 11. Gu, X.; Yang, R. First-Principles Prediction of Phononic Thermal Conductivity of Silicene: A Comparison with Graphene. J. Appl. Phys. 2015, 117, 025102. 27



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 28 of 49

12. Ward, A.; Broido, D.; Stewart, D. A.; Deinzer, G. Ab Initio Theory of the Lattice Thermal Conductivity in Diamond. Phys. Rev. B 2009, 80, 125203. 13. Broido, D.; Malorny, M.; Birner, G.; Mingo, N.; Stewart, D. Intrinsic Lattice Thermal Conductivity of Semiconductors from First Principles. Appl. Phys. Lett. 2007, 91, 231922. 14. Kong, B. D.; Paul, S.; Nardelli, M. B.; Kim, K. W. First-Principles Analysis of Lattice Thermal Conductivity in Monolayer and Bilayer Graphene. Phys. Rev. B 2009, 80, 033406. 15. Li, W.; Carrete, J.; Mingo, N. Thermal Conductivity and Phonon Linewidths of Monolayer MoS2 from First Principles. Appl. Phys. Lett. 2013, 103, 253103. 16. Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB , a Sparse Matrix-Based Implementation of the DFTB Method. J. Phys. Chem. A 2007, 111, 5678-5684. 17. Koskinen, P.; Mäkinen, V. Density-Functional Tight-Binding for Beginners. Comput. Mater. Sci. 2009, 47, 237-253. 18. Frauenheim, T.; Seifert, G.; Elsterner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. A Self-Consistent Charge Density-Functional Based Tight-Binding Method for Predictive Materials Simulations in Physics, Chemistry and Biology. Phys. Status Solidi B 2000, 217, 41-62. 19. Wang, X.; Guo, R.; Xu, D.; Chung, J.; Kaviany, M.; Huang, B. Anisotropic Lattice Thermal Conductivity and Suppressed Acoustic Phonons in MOF-74 from First Principles. J. Phys. Chem. C 2015, 119, 26000-26008. 20. Kuang, Y.; Lindsay, L.; Huang, B. Unusual Enhancement in Intrinsic Thermal Conductivity of Multilayer Graphene by Tensile Strains. Nano Lett. 2015, 15, 6121-6127. 21. Sieck, A.; Frauenheim, T.; Jackson, K. Shape Transition of Medium‐sized Neutral Silicon Clusters. Phys. Status Solidi B 2003, 240, 537-548. 22. Glassbrenner, C.; Slack, G. A. Thermal Conductivity of Silicon and Germanium from 3 K to the Melting Point. Phys. Rev. 1964, 134, A1058.

28


Page 29 of 49

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


23. Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260. 24. Grundkötter-Stock, B.; Bezugly, V.; Kunstmann, J.; Cuniberti, G.; Frauenheim, T.; Niehaus, T. A. SCC-DFTB Parametrization for Boron and Boranes. J. Chem. Theory Comput. 2012, 8, 1153-1163. 25. Kubař, T.; Bodrog, Z.; Gaus, M.; Köhler, C.; Aradi, B.; Frauenheim, T.; Elstner, M. Parametrization of the SCC-DFTB Method for Halogens. J. Chem. Theory Comput. 2013, 9, 2939-2949. 26. Elstner, M. The SCC-DFTB Method and its Application to Biological Systems. Theor. Chem. Acc. 2006, 116, 316-325. 27. Frenzel, J.; Oliveira, A. F.; Jardillier, N.; Heine, T.; Seifert, G. Semi-Relativistic, Self-Consistent Charge Slater–Koster Tables for Density-Functional Based Tight-Binding (DFTB) for Materials Science Simulations. http://www.dftb.org/parameters/download/matsci/matsci_0_3 (accessed on June 1, 2016). 28. Perdew, J. P.; Zunger, A. Self-Interaction Correction to Density-Functional Approximations for Many-Electron Systems. Phys. Rev. B 1981, 23, 5048. 29. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation made Simple. Phys. Rev. Lett. 1996, 77, 3865. 30. Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953. 31. Birch, F. Finite Elastic Strain of Cubic Crystals. Phys. Rev. 1947, 71, 809. 32. Cahangirov, S.; Topsakal, M.; Ciraci, S. Armchair Nanoribbons of Silicon and Germanium Honeycomb Structures. Phys. Rev. B 2010, 81, 195120. 33. Mehl, M. J. Pressure Dependence of the Elastic Moduli in Aluminum-Rich Al-Li Compounds. Phys. Rev. B 1993, 47, 2493.

29



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 30 of 49

34. Fugallo, G.; Lazzeri, M.; Paulatto, L.; Mauri, F. Ab Initio Variational Approach for Evaluating Lattice Thermal Conductivity. Phys. Rev. B 2013, 88, 045430. 35. Turney, J.; Landry, E.; McGaughey, A.; Amon, C. Predicting Phonon Properties and Thermal Conductivity from Anharmonic Lattice Dynamics Calculations and Molecular Dynamics Simulations. Phys. Rev. B 2009, 79, 064301. 36. Guo, R.; Wang, X.; Kuang, Y.; Huang, B. First-Principles Study of Anisotropic Thermoelectric Transport Properties of IV-VI Semiconductor Compounds SnSe and SnS. Phys. Rev. B 2015, 92, 115202.

37. Guo, R.; Wang, X.; Huang, B. Thermal Conductivity of Skutterudite CoSb3 from First Principles: Substitution and Nanoengineering Effects. Sci. Rep. 2015, 5. 38. O'Mara, W.; Herring, R. B.; Hunt, L. P. Handbook of Semiconductor Silicon Technology; Crest Publishing House, 2007. 39. Csonka, G. I.; Perdew, J. P.; Ruzsinszky, A.; Philipsen, P. H.; Lebègue, S.; Paier, J.; Vydrov, O. A.; Ángyán, J. G. Assessing the Performance of Recent Density Functionals for Bulk Solids. Phys. Rev. B 2009, 79, 155107. 40. Hopcroft, M. A.; Nix, W. D.; Kenny, T. W. What is the Young's Modulus of Silicon?. J. Microelectromech. S. 2010, 19, 229-238. 41. Haas, P.; Tran, F.; Blaha, P. Calculation of the Lattice Constant of Solids with Semilocal Functionals. Phys. Rev. B 2009, 79, 085104. 42. Park, J.; Yu, B. D.; Hong, S. Van Der Waals Density Functional Theory Study for Bulk Solids with BCC, FCC, and Diamond Structures. Curr. Appl. Phys. 2015, 15, 885-891. 43. Nilsson, G.; Nelin, G. Study of the Homology between Silicon and Germanium by Thermal-Neutron Spectrometry. Phys. Rev. B 1972, 6, 3777. 44. Dolling, G.; Cowley, R. The Thermodynamic and Optical Properties of Germanium, Silicon, Diamond and Gallium Arsenide. Proc. Phys. Soc. London 1966, 88, 463.

30


Page 31 of 49

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


45. Wei, S.; Chou, M. Phonon Dispersions of Silicon and Germanium from First-Principles Calculations. Phys. Rev. B 1994, 50, 2221. 46. Porter, L. J.; Justo, J. F.; Yip, S. The Importance of Grüneisen Parameters in Developing Interatomic Potentials. J. Appl. Phys. 1997, 82, 5378-5381. 47. Lalmi, B.; Oughaddou, H.; Enriquez, H.; Kara, A.; Vizzini, S.; Ealet, B.; Aufray, B. Epitaxial Growth of a Silicene Sheet. Appl. Phys. Lett. 2010, 97, 223109. 48. Tao, L.; Cinquanta, E.; Chiappe, D.; Grazianetti, C.; Fanciulli, M.; Dubey, M.; Molle, A.; Akinwande, D. Silicene Field-Effect Transistors Operating at Room Temperature. Nat. Nanotechnol. 2015, 10, 227-231. 49. Zhang, X.; Xie, H.; Hu, M.; Bao, H.; Yue, S.; Qin, G.; Su, G. Thermal Conductivity of Silicene Calculated using an Optimized Stillinger-Weber Potential. Phys. Rev. B 2014, 89, 054310. 50. Zhang, X.; Bao, H.; Hu, M. Bilateral Substrate Effect on the Thermal Conductivity of Two-Dimensional Silicon. Nanoscale 2015, 7, 6014-6022. 51. Kuang, Y. D.; Lindsay, L.; Shi, S. Q.; Zheng, G. Tensile Strains Give Rise to Strong Size Effects for Thermal Conductivities of Silicene, Germanene and Stanene. Nanoscale 2016, 8, 3760-3767. 52. Peng, B.; Zhang, H.; Shao, H.; Xu, Y.; Ni, G.; Zhang, R.; Zhu, H. Phonon Transport Properties of Two-Dimensional Group-IV Materials from Ab Initio Calculations. Phys. Rev. B 2016, 94, 245420. 53. Cahangirov, S.; Topsakal, M.; Aktürk, E.; Şahin, H.; Ciraci, S. Two-and One-Dimensional Honeycomb Structures of Silicon and Germanium. Phys. Rev. Lett. 2009, 102, 236804. 54. Kuang, Y.; Lindsay, L.; Shi, S.; Wang, X.; Huang, B. Thermal Conductivity of Graphene Mediated by Strain and Size. Int. J. Heat Mass Transfer 2016, 101, 772-778. 55. Xie, H.; Ouyang, T.; Germaneau, É; Qin, G.; Hu, M.; Bao, H. Large Tunability of Lattice Thermal Conductivity of Monolayer Silicene Via Mechanical Strain. Phys. Rev. B 2016, 93, 075404. 31



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 32 of 49

56. Esfarjani, K.; Stokes, H. T. Method to Extract Anharmonic Force Constants from First Principles Calculations. Phys. Rev. B 2008, 77, 144112. 57. Liu, Z.; Wu, X.; Luo, T. The Impact of Hydrogenation on the Thermal Transport of Silicene. 2D Mater. 2017, 4, 025002.

58. Carrete, J.; Li, W.; Lindsay, L.; Broido, D. A.; Gallego, L. J.; Mingo, N. Physically Founded Phonon Dispersions of Few-Layer Materials and the Case of Borophene. Mater. Res. Lett. 2016, , 1-8. 59. Ma, J.; Chen, Y.; Han, Z.; Li, W. Strong Anisotropic Thermal Conductivity of Monolayer WTe2. 2D Mater. 2016, 3, 045010. 60. Topsakal, M.; Ciraci, S. Elastic and Plastic Deformation of Graphene, Silicene, and Boron Nitride Honeycomb Nanoribbons Under Uniaxial Tension: A First-Principles Density-Functional Theory Study. Phys. Rev. B 2010, 81, 024107.

32


Page 33 of 49

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


Figure 1. Calculated phonon dispersion relations of silicon from DFT/PBE method, DFTB method with pbc and matsci sets along with the experiment values43,44. Phonon DOS from the DFT/PBE method and DFTB method with pbc and SiTherm sets are also shown on the right panel.

Figure 2. Calculated phonon dispersion relations of silicon from MD with SW and Tersoff potentials comparing with those from DFTB/SiTherm and experiment values43,44.

33



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 34 of 49

Figure 3. Calculated mode Grüneisen parameters of silicon from the DFT method and the DFTB method with SiTherm set along with experiment values5.

Figure 4. Calculated thermal conductivities of silicon together with experiment results22 from 200 K to 700 K. Stars are DFT results from Ref. 13.

34


Page 35 of 49

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


Figure 5. Calculated phonon mode relaxation times in silicon crystal.

Figure 6. Optimized structures of 9ASiNR and 7ZSiNR from the DFT/LDA [(a) and (c)] and DFTB/SiTherm [(b) and (d)] methods. Edge reconstructions32,53 are well reproduced.

35



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 36 of 49

Figure 7. Calculated phonon dispersion relations and DOS of silicene. The inset is the logarithmic plot of the dispersion relation near the Γ point.

Figure 8. Calculated lattice thermal conductivities of silicene from the DFT/LDA and DFTB/SiTherm methods from 200 K to 800 K. 36


Page 37 of 49

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


Figure 9. Calculated phonon relaxation times in silicene from the DFT/LDA and DFTB/SiTherm methods.

37



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 38 of 49

Table 1. Comparison of lattice constants, cohesive energy and elastic moduli of silicon predicted by the DFT and DFTB methods and experimental results38-40.

Lattice

DFTB/SiTherm DFTB/pbc

DFTB/matsci27 DFT/LDA DFT/PBE Experiment

5.421

5.460

5.147

5.403

5.469

5.430738

4.12

4.83

7.69

5.35

4.60

4.69339

157.5

153.2

285.5

154.9

152.2

130-18840

60.6

61.8

102.0

63.0

62.7

51-8040

constant (Å) Cohesive energy (eV) Young’s modulus (GPa) Shear Modulus (GPa)

38


Page 39 of 49

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


Table 2. Computed structural parameters of silicene, in Å. Other DFT data are also listed in parentheses. Silicene

Lattice constant

Buckling height

(Å)

(Å)

DFT/LDA

3.826 (3.82551)

0.438 (0.43951)

DFT/PBE

3.868 (3.8752)

0.450 (0.4452)

DFTB/pbc

3.814

0.617

DFTB/matsci

3.820

0

DFTB/SiTherm

3.816

0.448

39



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 40 of 49

TOC Graphic

40


Page 41 of 49

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


Figure 1. Calculated phonon dispersion relations of silicon from DFT/PBE method, DFTB method with pbc and matsci sets along with the experiment values41,42. Phonon DOS from the DFT/PBE method and DFTB method with pbc and SiTherm sets are also shown on the right panel. 59x42mm (300 x 300 DPI)



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

Figure 2. Calculated phonon dispersion relations of silicon from MD with SW and Tersoff potentials comparing with DFTB method with SiTherm set and experiment values41,42. 65x52mm (300 x 300 DPI)


Page 42 of 49

Page 43 of 49

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


Figure 3. Calculated mode Grüneisen parameters of silicon from DFT method and DFTB method with SiTherm set along with experiment values5. 63x48mm (300 x 300 DPI)



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

Figure 4. Calculated thermal conductivities of silicon with experiment results22 from 200 K to 700 K. Stars are DFT results from Ref. 13. 62x47mm (300 x 300 DPI)


Page 44 of 49

Page 45 of 49

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


Figure 5. Calculated phonon mode relaxation times in silicon crystal. 59x43mm (300 x 300 DPI)



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

Figure 6. Optimized structures of 9ASiNR and 7ZSiNR from the DFT/LDA [(a) and (c)] and DFTB/SiTherm [(b) and (d)] methods. Edge reconstructions50,51 are well reproduced. 84x86mm (300 x 300 DPI)


Page 46 of 49

Page 47 of 49

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


Figure 7. Calculated phonon dispersion relations and DOS of silicene. The inset is the logarithmic plot of the dispersion relation near the Γ point. 60x38mm (300 x 300 DPI)



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

Figure 8. Calculated lattice thermal conductivities of silicene from DFT/LDA and DFTB/SiTherm methods from 200 K to 800 K. 62x47mm (300 x 300 DPI)


Page 48 of 49

Page 49 of 49

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


Figure 9. Calculated phonon relaxation times in silicene from the DFT/LDA and DFTB/SiTherm methods. 60x44mm (300 x 300 DPI)


Parametrization of Density Functional Tight-Binding Method for

Recommend Documents