Home

Search

Collections

Journals

About

Contact us

My IOPscience

Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 J. Phys.: Condens. Matter 26 225402 (http://iopscience.iop.org/0953-8984/26/22/225402) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 132.174.255.116 This content was downloaded on 17/06/2014 at 12:37

Please note that terms and conditions apply.

Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 26 (2014) 225402 (12pp)

doi:10.1088/0953-8984/26/22/225402

Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations T Tadano1, Y Gohda1,2 and S Tsuneyuki1,3 1

  Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan   Department of Materials Science and Engineering, Tokyo Institute of Technolgy, Yokohama 226-8502, Japan 3   Institute for Solid State Physics, The University of Tokyo, Kashiwa 277-8581, Japan 2

E-mail: [email protected] Received 1 January 2014, revised 24 March 2014 Accepted for publication 28 March 2014 Published 14 May 2014 Abstract

A systematic method to calculate anharmonic force constants of crystals is presented. The method employs the direct-method approach, where anharmonic force constants are extracted from the trajectory of first-principles molecular dynamics simulations at high temperature. The method is applied to Si where accurate cubic and quartic force constants are obtained. We observe that higher-order correction is crucial to obtain accurate force constants from the trajectory with large atomic displacements. The calculated harmonic and anharmonic force constants are, then, combined with the Boltzmann transport equation (BTE) and non-equilibrium molecular dynamics (NEMD) methods in calculating the thermal conductivity. The BTE approach successfully predicts the lattice thermal conductivity of bulk Si, whereas NEMD shows considerable underestimates. To evaluate the linear extrapolation method employed in NEMD to estimate bulk values, we analyze the size dependence in NEMD based on BTE calculations. We observe strong nonlinearity in the size dependence of NEMD in Si, which can be ascribed to acoustic phonons having long mean-free-paths and carrying considerable heat. Subsequently, we also apply the whole method to a thermoelectric material Mg2Si and demonstrate the reliability of the NEMD method for systems with low thermal conductivities. Keywords: anharmonic force constant, lattice thermal conductivity, non-equilibrium molecular dynamics, Boltzmann transport equation, first-principles calculation S Online supplementary data available from stacks.iop.org/J.PhysCM/26/225402/mmedia (Some figures may appear in colour only in the online journal)

conductivity and increasing the power factor σS2. Since the heat is mainly carried by lattice vibrations in typical thermoelectric materials, ideal thermoelectric materials follow the phonon-glass electron-crystal (PGEC) concept proposed by Slack [1]. Since the theoretical predictions by Hicks and Dresselhaus [2, 3], nanostructures have been successfully introduced to reduce κℓ. For example, Bi2Te3/Sb2Te3 superlattice was reported to show ZT ∼ 2.4 [4], and recent experimental studies showed that silicon nanowires with rough surfaces have significantly higher ZT compared to the bulk counterpart [5, 6]. In order to understand the mechanism of reduction in κℓ in thermoelectric materials and to further improve the

1. Introduction Thermal conductivity has received much attention in recent years because of its importance in micro-/nanostructured semiconducting devices. Thermoelectric materials are one of important targets, which are expected to be a next-generation device for energy conversion. The efficiency of a thermoelectric material is given by the dimensionless figure of merit ZT = σ S2T/(κc+κℓ), where σ is the electrical conductivity, S is the Seebeck coefficient, T is the absolute temperature, and κc (ℓ) is the thermal conductivity of electrons (phonons), respectively. ZT can be enhanced by reducing the thermal 0953-8984/14/225402+12$33.00

1

© 2014 IOP Publishing Ltd  Printed in the UK

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

that is considered to be a source of uncertainty in κℓ [20]. To emphasize the validity of our computational approach and elucidate the size effect in NEMD, we also apply the method to Mg2Si, which is a thermoelectric material having much lower κℓ.

efficiency, reliable theoretical/computational approaches for the lattice thermal conductivity are highly demanded. There are several methods to calculate the lattice thermal conductivity such as the Boltzmann transport equation (BTE), equilibrium molecular dynamics (MD) with the Green–Kubo formula (EMD-GK) and non-equilibrium MD (NEMD). Each of them has been widely employed in combination with empirical model potentials [7–11]. In BTE, it is easier to analyze mode-dependent transport properties of phonons than in MD methods. However, forth- and higher-order anharmonicity are usually neglected in BTE, so that the scattering probability of phonons can be underestimated at high temperatures. MD methods, on the other hand, consider full anharmonicities and are applicable to complex structures with surfaces and interfaces if any reliable model potential is available [9, 10]. Since well-established empirical potential is limited, however, it has been difficult to estimate and analyze phonon transport in complex thermoelectric materials. To overcome the limitation, first-principles approaches have recently been employed for both BTE [12–14], and MD methods [15]. For example, Broido et al [12] estimated the lattice thermal conductivity of specific bulk materials using BTE with harmonic and third-order anharmonic force constants obtained by ab initio calculations based on density-functional theory (DFT), and the results showed a good agreement with experimental values. Although it is desirable to use a firstprinciples method instead of empirical model potentials in MD simulations and there is an example [15], a straightforward application is formidable in most cases due to its high computational costs attributed to the large system size and the long simulation time required for well-converged results. One way to handle this difficulty is to employ first-principles calculations in an indirect manner, in which the parameters of model potentials are determined from the results of DFT calculations and large-scale MD simulations are performed with optimized model potentials as in references [16, 17]. In the present paper, we construct versatile model potentials from first principles that enable us to predict the lattice thermal conductivity of solids either by BTE or MD. The model potential is represented as a Taylor expansion with respect to atomic displacements, where harmonic and anharmonic force constants are the parameters to be determined. In this study, we calculate anharmonic force constants by using first-principles molecular dynamics (FPMD) instead of the conventional technique where the multiple atoms in the supercell are displaced manually [18, 19]. In FPMD, probable patterns of ionic displacements are automatically generated and the sampling space is easily controlled through the temperature. In section 2, we introduce the interatomic model potential for thermal conductivity predictions and the detailed procedure to determine parameters from first principles. To demonstrate the validity of our approach, we apply the methodology to Si, a typical covalent semiconductor. Here, we investigate the dependencies of anharmonic force constants on temperature and computational conditions of FPMD. In section 3, we apply harmonic and anharmonic force constants to thermal conductivity calculations. Here, we employ both BTE and NEMD and compare these results. In addition, we also investigate the size effect in NEMD simulations,

2.  Interatomic model potential and parameter fitting 2.1.  Anharmonic lattice model

In solids, it is usual that the atoms oscillate about their equilibrium positions. As long as the displacements are small compared to the interatomic distances, it is justified to approximate the interatomic potential as a Taylor expansion with respect to atomic displacements from equilibrium positions at zero temperature: V − V0 =



∑

N =2

1 N!

∑

∑

m1α1, ⋯ , mN αN μ1, ⋯ , μN

× u mμ11α1 ⋯

μ1, ⋯ , μN Φm 1α1, ⋯ , mN αN

u mμNN αN .

(1)

Here, mi is a translation vector of the primitive lattice, αi is the atom index in the unit cell at mi, and μi specifies xyz components, respectively. u mμα is the displacement of atom α in the unit cell m along the μ direction and the coupling constants Φ's are interatomic force constants (IFCs) which are defined by ∂N V μ1, ⋯ , μN Φm . 1α1, ⋯ , mN αN = (2) ∂ u mμ11α1 ⋯∂ u mμNN αN u mμ1α =⋯= 0 1 1

The model defined by equation (1), which we call anharmonic lattice model (ALM), has two control parameters. One is the maximum order Nmax which specifies where to truncate the series, and the other is the cutoff radius rN which restricts the interacting atoms (m1α1, ···, mNαN) in the summation. IFCs are parameters of the ALM which depend on atomic species and structures. The number of parameters increases rapidly with rN and Nmax. In practical calculations, relations listed below can be used to reduce the number of independent parameters. Permutation.  From the definition (equation (2)), force con-

stants are invariant under an exchange of any two sets of index pairs. That is, μ1, μ2 , μ3, ⋯ μ1, μ3, μ2 , ⋯ Φm (3) 1α1, m2 α2, m3 α3, ⋯ = Φ m1α1, m3 α3, m2 α2, ⋯ = ⋯ ,

which can reduce the number of parameters roughly by 1/N! for Nth-order IFCs. Periodicity.  Because the IFCs depend on interatomic dis-

tances, they are invariant under a translation in units of vector m, which can be formulated as μ1, μ2 , μ3, ⋯ μ1, μ2 , μ3, ⋯ Φm (4) 1α1, m2 α2, m3 α3, ⋯ = Φ 0α1, m2 − m1α2, m3 − m1α3, ⋯ = ⋯ .

Therefore, even if atom α1 is restricted to be within a specific unit cell, say 0, all independent parameters can be completely determined. Crystal symmetry.  A symmetry operation of a structure relates an equilibrium atomic position xmα to another one x m′α′ 2

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

where the first term is the Hellmann–Feynman force acting on the atom mα at ith configuration and the second term is that of the model force which can be obtained by partially differentiating the ALM by umα. Here, M is the number of displaced configurations used as the reference. The minimization of χ2 results in a linear least-square fitting with the equality constraints (equation (6)) which can be directly solved without iteration. The number of independent data set has to be equal to or greater than the number of the independent parameters; otherwise the problem is underdetermined and the parameters cannot be determined uniquely. The accuracy of the model for the reference system can be estimated with the relative error σ given by

through rotation and translation. If the system has Ns such symmetry operations, there are Ns equations which relates an IFC to another one by 

∑

ν1, ν2, ⋯ , νN

, νN Φνm1, ′να2,′⋯ , m ′α ′, ⋯ , m 1 1

2

2

N

′

αN ′

Oμ1ν1Oμ2ν2 ⋯ OμN νN

μ1, μ2 , ⋯ , μN = Φm 1α1, m2 α2, ⋯ , mN αN ,

(5)

where O is the rotation matrix of the symmetry operation. In addition to the relations listed above, parameters need to satisfy constraint relations for the translational/rotational symmetry of ALMs. Especially, in three dimensional periodic systems, the translational invariance relation is necessary to obtain correct phonon dispersion relations, which is given by

χ2 μ1, μ2 , ⋯ , μN ∑ Φ m σ = , 1α1, m2 α2, ⋯ , mN αN = 0 , ∀ ( m2 α2, ⋯ , mN αN ; μ1 , μ 2 , ⋯ , μN ) . M (8) m1α1 ∑i ∑mα Fmα ( ti ) 2 (6)

At finite temperatures, the center of atomic oscillation can be slightly changed due to thermal expansion. We find that, with the invariance relations, ALMs can simulate the effects of the thermal expansion even though they are defined at the absolute zero temperature (equation (1)). The details of this property will be discussed elsewhere.

where the residual sum of squares is normalized by the square sum of Hellmann–Feynman forces. To generate data sets for the fitting, it is usual that atoms in the supercell are manually displaced and Hellmann–Feynman forces are calculated for each configuration. In the case of harmonic force constants, the magnitude of the displacement Δu has to be small so that anharmonic contributions can be neglected, which is typically of the order of 0.01 Å. In practical calculations, multiple data sets with different magnitude of displacements are usually used to reduce the effects of anharmonicity and numerical errors. For harmonic IFCs, this procedure has been employed successfully. For anharmonic cases, much attention needs to be paid to extract reliable IFCs. Typically, the potential energy surface is almost harmonic, and for small displacements, say 0.01 Å, anharmonic corrections to the forces can be as small as the errors in Hellmann–Feynman forces. In such cases, the accuracy of DFT calculations needs to be improved by employing finer k-point grids, smaller convergence thresholds in the self-consistent loop, and so forth. Alternatively, the errors in anharmonic IFCs can be reduced by employing large displacements. The larger the displacement, the more significant the anharmonic forces, and so the sampling can be carried out more accurately. In order to obtain data sets for Nth-order IFCs, N − 1 atoms at most need to be displaced along different directions simultaneously. The number of such patterns increases rapidly when N becomes large, which must be tedious to handle manually. In order to extract anharmonic IFCs in a simpler manner, we alternatively employed first-principles molecular dynamics (FPMD) in this report. In FPMD simulations, all atoms in the supercell displace from equilibrium positions at zero temperature. Therefore, once the displacement-force data are obtained, it can be used to extract IFCs of an arbitrary order. Also, the probable patterns of displacements are easily obtained and the magnitude of displacements can be controlled through the temperature. Furthermore, by virtue of the properties defined by equations (4) and (5), it is possible to sample several independent data sets from each configuration of FPMD simulation, which can help reduce computational costs associated with the data sampling. In the present study,

2.2.  First-principles calculations of force constants

There are two different methods to calculate IFCs from firstprinciples: density functional perturbation theory (DFPT) and the direct method. In DFPT, first-order changes of the wave function are calculated self-consistently [21, 22]. Once the firstorder change is obtained, energy perturbations of the second and third orders can be calculated, which correspond to harmonic and cubic IFCs, respectively. With the first-order DFPT, harmonic and cubic force constants of arbitrary k-points can be calculated efficiently with a primitive cell. However, it has a difficulty in calculating higher-order energy perturbations where the higher-order changes of the wave function are required. In the direct method, IFCs are estimated via Hellmann–Feynman forces acting on atoms in displaced configurations [23]. If distant interactions are essential, supercell calculations will be necessarily for reliable IFCs. One advantage of the direct method is that anharmonic IFCs of fourth and higher orders can be easily extracted in a systematic manner, which is formidable in DFPT. Also, the importance of distant interactions can be easily examined through the magnitude of real-space harmonic and anharmonic IFCs. Therefore, we choose the direct method approach to extract all parameters in this study. Upon calculating harmonic IFCs in the direct method, it is sufficient to displace just one atom from its equilibrium position, and the forces acting on the atoms are used. The same strategy can be used to calculate anharmonic force constants where multiple atoms in the supercell need to be displaced simultaneously. For both harmonic and anharmonic cases, we estimate IFCs numerically via minimization of the residual sum of squares of the forces: M

2 ∼ χ 2 = ∑∑ Fmα ( ti ) − Fmα ( ti ) , (7) i

mα

3

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

we assume that effects of thermal excitations of electrons on IFCs are negligible for semiconductors. 2.3.  Computational conditions

To see the validity of the strategy introduced above, bulk silicon was chosen as a typical semiconductor. In order to consider distant interaction pairs, a 3  ×  3  ×  3 conventional cell with 216 atoms was prepared. For the sampling of displacement-force data sets, first-principles calculations were carried out using QUANTUM-ESPRESSO package [24], which uses the plane-wave-pseudopotential method based on DFT. A norm-conserving pseudopotential was used and the cutoff energy of 32  Ry was employed for wave functions. For the exchange-correlation functional, we employed local density approximation (LDA) with Perdew–Zunger parametrization [25, 26]. Unless otherwise mentioned, 3 × 3 × 3 MonkhorstPack k-point grid was used for the Brillouin zone integration [27]. Before sampling the data, the cell parameter was optimized and we have obtained the lattice constant of 5.40  Å, comparable with the experimental value 5.43 Å. Upon generating the reference data set for harmonic IFCs, an atom in the supercell was slightly displaced along the [1 0 0] direction. To reduce numerical errors, the amount of displacements were chosen to be ±Δu, ±2Δu, and ±3Δu, where Δu  =  0.008  Å, and the harmonic IFCs were extracted from these six data sets. For anharmonic IFCs, 2000 steps of FPMD were performed, where first 500 data were discarded because they retained the thermalizing nature. The temperature was controlled with Berendsen thermostat [28], and the time step of 0.96 fs was employed. To see the dependence on the control temperature and the fineness of the k-point grid, we prepared three different references. One of them was generated by a FPMD simulation with 3 × 3 × 3 k points at 1000 K, another is a gamma point calculation at the same temperature, and the other is a gamma point calculation at 500 K, where they will be referred to as FPMD1, FPMD2, and FPMD3, respectively. Upon extracting anharmonic IFCs, harmonic IFCs were fixed. In both harmonic and anharmonic calculations, numerical fittings were performed by the QR-decomposition method implemented in the DGGLSE routine of LAPACK.

Figure 1. Phonon dispersion relations of bulk Si calculated from harmonic IFCs. The dotted (black), dashed (blue), and thin solid (green) lines are results of second-, fifth-, and ninth-NN models, respectively. The results of a DFPT calculation (thick solid line) and experimental measurements by Nilsson and Nelin [29] (triangle) and by Kulda et al [30] (circle) are also shown for comparison.

radius. Aouissi et al [31] investigated the effect of the interaction cutoff on the phonon dispersion and reported ωLA(L) of ninth-NN model as 364 cm−1, which is 363.9 cm−1 in this study, and that of 25th-NN shell as 374 cm−1, which agrees better with the experimental result 378  cm−1. Thus, if more accurate description of the phonon dispersion is required, more distant interactions should be included. Compared to the result of the DFPT calculation and the experimental result, the accuracy of ninth-NN interaction model is satisfactory (figure 1). Therefore, the ninth-NN interaction model is employed throughout the study. 2.5.  Anharmonic IFCs from FPMD

In this work, anharmonic terms are considered up to sixth order within the first-NN shell. The numbers of independent IFCs are 5, 14, 22, and 46 for third to sixth orders, respectively. In table 1, calculated cubic and quartic IFCs are summarized together with the results of the previous work [18]. In the table and the following, the index m in anharmonic IFCs are omitted for simplicity because all m's are identical in this study. The subscripts 1 and 2 in IFCs specify the atoms located at (0, 0, 0) and a (1/4, 1/4, 1/4) respectively, where a is the lattice constant. Cubic force constants extracted from different FPMD trajectories agree well with each other. Also, the agreement with the manual displacement method is quite satisfactory, in which atom 1 and atom 2 were simultaneously displaced along 〈1 1 1〉 directions with the amount of |Δu| = 0.0648 Å and the fourthorder ALM is employed as the fitting function. The error (equation (8)) in the fitting to the manually displaced data set is 0.0127. In the FPMD methods, the atomic displacements are far larger than those in the manual displacement, therefore we employ the sixth-order ALM for the numerical fitting. To compare the magnitude of displacements, the root-meansquare displacement (RMSD) of atoms is calculated from the configurations of FPMD as shown in table 1. Even though the

2.4.  Extraction of harmonic IFCs

For harmonic IFCs, interaction pairs up to ninth nearest neighbors (NNs) are considered, resulting in 38 independent parameters. The fitting error σ is 0.0173 that is reasonably small. In figure 1, the calculated phonon dispersion relation is shown together with that from a DFPT calculation with 4 × 4 × 4 k points and experimental measurements [29, 30]. In order to investigate the effect of the distant interactions, we also examined the cases where interaction pairs were restricted up to second- and fifth-NN shells with the fitting errors of 0.0383 and 0.0366, respectively. The corresponding phonon dispersion relations are also shown in figure 1. From the results of second- and fifth-NN models, a minor underestimation in the frequency of the longitudinal acoustic (LA) mode in the vicinity of the L point can be ascribed to the introduction of a cutoff 4

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402 −3

−4

Table 1.  Cubic (unit: eV Å ) and quartic (unit: eV Å ) force constants extracted from three different FPMD trajectories: FPMD1 (3 × 3 × 3 k-points, T = 1000 K), FPMD2 (Γ point, T = 1000 K) and FPMD3 (Γ point, T = 500 K). Results of the conventional displacement method and previous work are also shown for comparison. The 95% confidence intervals are shown for the results of FPMD. The index m of IFCs is constant here, hence it is omitted.

Sampling method RMSD (Å) x , y, z Φ1,1,1 x, x, x Φ1,1,2 x, x, y Φ1,1,2 x , y, x Φ1,1,2 x , y, z Φ1,1,2 x, x, x, x Φ1,1,1,1 x, x, y, y Φ1,1,1,1 x, x, x, x Φ1,1,1,2 x, x, x, y Φ1,1,1,2 x, x, y, x Φ1,1,1,2 x, x, y, y Φ1,1,1,2 x, x, y, z Φ1,1,1,2 x , y, z , x Φ1,1,1,2 x, x, x, x Φ1,1,2,2 x, x, x, y Φ1,1,2,2 x, x, y, y Φ1,1,2,2 x, x, y, z Φ1,1,2,2 x , y, x , y Φ1,1,2,2 x , y, x , z Φ1,1,2,2

a

Present work FPMD1 0.165

FPMD2 0.168

FPMD3 0.113

Displacement 0.0374a

Previous work [18] Displacement 0.008, 0.016a

31.69(4) −3.10(1) −6.07(1) −6.07(1) −7.92(1) −59.17 (43) 39.42 (25) 14.79 (11) −4.08(6) −4.08(6) −9.85(6) −16.27 (9) −16.27 (9) −14.79 (11) 4.08 (6) 9.85 (6) 16.27 (9) 9.85 (6) 16.27 (9)

31.64 (5) −3.13(2) −6.07(1) −6.07(1) −7.91(1) −60.38 (37) 38.69 (28) 15.09 (9) −3.81(5) −3.81(5) −9.67(7) −16.08 (9) −16.08 (9) −15.09 (9) 3.81 (5) 9.67 (7) 16.08 (9) 9.67 (7) 16.08 (9)

31.88 (3) −3.07(1) −6.09(1) −6.09(1) −7.97(1) −56.23 (45) 40.55 (32) 14.06 (11) −3.95(6) −3.95(6) −10.14 (8) −16.63 (9) −16.63 (9) −14.06 (11) 3.95 (6) 10.14(8) 16.63(9) 10.14(8) 16.63(9)

32.81 −2.94 −6.02 −6.02 −8.20 −48.63 45.32 12.16 −3.65 −3.65 −11.33 −17.69 −17.69 −12.16 3.65 11.33 17.69 11.33 17.69

32.26 −2.74 −6.29 −6.29 −8.06 −42.99 40.15 10.75 −2.99 −2.99 −10.04 −14.62 −14.62 −10.75 2.99 10.04 14.62 10.04 14.62

The magnitude of displacements taken along xyz directions.

the SW potential: one with the original parameter set, and the other with the optimized parameters which is based on LDA calculations [17]. The original SW potential tends to overestimate forces, which is partly improved in the optimized SW potential. Compared to these SW potentials, the accuracy of the fourth-order ALM is much better. Thus, reliable description of dynamical and transport properties of phonons is expected by first-principles ALMs. Fifth- and sixth-order force constants are also estimated at the same time from FPMD results at high temperatures. We observe that their contributions to the total potential energy is an order of magnitude smaller than those of cubic and quartic terms even at 1000 K, which makes it difficult to determine accurate higher-order IFCs. Even though fifth- and sixth-order terms have some uncertainties, they can be employed in MD simulations for the improved stability of ALMs at high temperatures. The fitting errors σ against 1500 MD trajectories are 0.0541, 0.0579 and 0.0378 for FPMD1, FPMD2 and FPMD3, respectively. The errors are not reduced by seventh-, eighth-, and higher-order contributions because the number of adjustable parameter doesn't increase much due to the constraints (equation (6)). Therefore, the residual errors are likely to be ascribed to more distant interactions neglected here.

displacement is larger in 1000  K than in 500  K, it is found that the cubic force constants are not so much sensitive to the temperature in FPMD, i.e. differences in cubic IFCs are about 2% at most, provided that the sixth-order ALM is employed for the fitting function. On the other hand, when fifth- and sixth-order terms are neglected upon fitting, the anharmonic IFCs become more sensitive to the temperature in FPMD. For x , y, z example, when the fourth-order ALM is employed, Φ1,1,1 is −3 −3 estimated to be 29.24 eV Å and 30.96 eV Å for FPMD2 (1000 K) and FPMD3 (500 K), respectively. The overall difference is increased to about 5% and the extracted values are significantly underestimated compared to the result of the manual displacement. This implies that the neglected higherorder (fifth- and sixth-order) contributions are renormalized into the lower-order (third- and fourth-order) anharmonic terms affecting their accuracy. Therefore, the inclusion of higher order anharmonic terms is crucial to extract accurate lower-order anharmonic IFCs from FPMD configurations at high temperatures. In figure 2, we compare Hellmann–Feynman forces and corresponding model forces in various configurations obtained from FPMD at 1000  K. Here, we show the model force in ALM without anharmonicity (ALM2) and that with anharmonicity up to quartic order (ALM4). As clearly seen from the figure, ALM2 significantly deviates from the FPMD data due to absence of anharmonicity. With cubic and quartic terms, the accuracy of ALM improves greatly and the the Hellmann– Feynman forces are well reproduced in a wide force range. To compare the accuracy with an empirical model potential, we also show atomic forces of the Stillinger–Weber (SW) potential [32]. In the figure, we show two different results for

3.  Lattice thermal conductivity from first principles 3.1. Overview

BTE, EMD-GK, and NEMD methods have been extensively employed to calculate the lattice thermal-conductivity of solids, such as mantle minerals [15, 33], thermoelectric 5

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

understand the size effects in NEMD and assess the quantitative reliability of NEMD method, thermal conductivity predictions based on first-principles should be informative. In the following, we first estimate κℓ of bulk Si both with BTE-RTA and NEMD. Subsequently, we also show computational results on a thermoelectric material Mg2Si. For both systems, the agreement between BTE and NEMD is evaluated based on first-principles force constants. 3.2.  BTE simulation

In this study, BTE within the single-mode relaxation time approximation is employed, in which the lattice thermal conductivity is given by 1 κ ℓμ, ν = ∑ c ( q, n ) v μ ( q, n ) v ν ( q, n ) τ ( q, n ) . (9) Ω Nq q, n

Here, c (q, n), v(q, n), and τ (q, n) are the constant-volume specific heat, the group velocity, and the relaxation time of phonons with wave vector q and polarization n, respectively. Ω is the volume of the unit cell and Nq is the number of q points. In defect-free bulk samples, three phonon interactions dominate the phonon scattering. Therefore, we approximate the relaxation time as τ−1(q, n)  ≈  2Γqn(ωqn), where Γqn(ω) is the imaginary part of the phonon selfenergy considering the lowest-order three-phonon process. Upon calculating Γqn(ω), we employ the tetrahedron method instead of broadening approaches [40]. For details of Γqn(ω) employed in this work and its integration, see Appendix A. Before moving on to the comparison with experimental results, we have investigated the convergence of lattice thermal-conductivity calculated by equation (9) with respect to the q-grid density and the cutoff radius of cubic IFCs. In the previous section, we extracted cubic IFCs restricted within the first-NN shell (1NN model). Here, in order to examine the importance of distant cubic interactions, we extended the cutoff radius so that cubic interactions were considered up to second- and third-NN shells (2NN and 3NN models). The number of cubic IFCs were found to be 36 and 95 for the 2NN and 3NN models, which were five in the 1NN model, and the anharmonic IFCs were extracted from the trajectory of FPMD, specifically that of FPMD2, in exactly the same manner as explained in section 2. Lattice thermal-conductivities of these three models calculated by equation (9) are compared in figure 3 with different number of q points. Firstly, we can see that the 1NN model underestimates κℓ of the 2NN and 3NN models by about 10% . Because the 1NN model has only five cubic IFCs and the translational invariance (equation (6)) gives rise to 2 constraint relations between them, the degree of freedom reduces to 3 that is found to be too small to optimize the cubic interaction in Si. When we consider cubic interactions beyond the first-NN shell, the number of parameters increase sufficiently to capture the anharmonic interaction. Since the anharmonic interaction is short-ranged, we observed no significant difference between the 2NN and 3NN models. Therefore, it is sufficient to consider cubic interactions up to the second-NN shell for

Figure 2. Comparison of first-principles forces and model forces in various configurations at 1000 K. The harmonic model (ALM2, circle) significantly deviates from the perfect match dashed-line due to lack of anharmonicity, which is greatly improved by the fourthorder model (ALM4, triangle). ALM4 reproduces first-principles forces better than the Stillinger–Weber potential with original (cross) and optimized (square) parameters.

materials [13, 14, 16], and nanowires [9, 10]. Among them, bulk Si has been intensively studied as a typical semiconductor [7, 11, 12, 34, 35]. For calculating κℓ of bulk Si, empirical model potentials, such as SW and Tersoff potentials [36], have been repeatedly employed due to their computational feasibility. According to the previous BTE and MD simulations, however, the accuracy of these empirical model potentials is unsatisfactory in predicting κℓ. For example, a NEMD simulation with SW potential yielded κℓ = 119 ± 40 Wm−1K−1 at 500 K [11], which overestimates an experimental value of 80  Wm−1K−1 [37]. Another study by BTE with SW potential showed considerably overestimated κℓ [7], i.e. κ ℓSW ≈ 900  Wm−1K−1 at 300 K that is about six times as large as an experimental value of 156 Wm−1K−1 [38]. These results suggest inaccurate vibrational properties of empirical potentials, which now can be improved by first-principles approaches. Broido et al [12] solved BTE with accurate harmonic and anharmonic IFCs obtained from DFT calculations and achieved great agreements with experimental results. More recently, Esfarjani et al [39] reported a detailed analysis on the phonon transport in bulk Si, where they employed BTE with the relaxation time approximation (RTA). Even if exactly the same interatomic potential is employed, considerable disagreements between BTE and NEMD have been observed for bulk Si. For the wide temperature range, bulk thermal conductivity of NEMD tends to underpredict that of BTE. Sellan et al [20] analyzed the size effect in NEMD simulations with the SW potential and suggested that the source of the underestimation can be the linear extrapolation approach employed in NEMD [11]. In order to further 6

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

Figure 3. Lattice thermal conductivity of bulk Si at 300 K calculated by BTE-RTA as a function of the number of q points. Cubic interactions are restricted within first-, second-, and thirdNN shells in 1NN, 2NN, and 3NN models, respectively. The linear extrapolation to the infinite q-grid density is shown by dashed lines. In addition, an experimental thermal conductivity of isotopically pure sample is indicated by the horizontal line.

Figure 4. Temperature dependence of the lattice thermal conductivity of bulk Si calculated by BTE-RTA (solid lines) and NEMD (circle with error bars) simulations together with the experimental values for isotopically enriched [38] (open circle) and naturally occurring isotope concentrations [37] (triangle). Previous computational results by BTE-RTA [39] (blue dotted line) and BTE [41] (green dashed line) are also shown for comparison.

predicting κℓ of Si. Secondly, we can see that κℓ increases as increasing the q-grid density. A denser q grid enables longwavelength phonons having large τ (q, n) to take part in heat transfer, which consequently increases the lattice thermalconductivity. In order to estimate κℓ in the infinite q density limit, we employed the analytic formula derived by Esfarjani and co-authors [39]. They showed that, assuming that the density of states and the relaxation time of long-wavelength acoustic phonons can be well described by the Debye model and the Klemens' formula respectively, the size-dependent thermal conductivity is given by κ ℓ ( Nqμ ) = κ ℓ∞− a / Nqμ with a proportionality constant a. As clearly seen in the figure, the size dependence follows the linear relationship in Nqμ ≳ 10. Extrapolation to infinite q-grid density ( Nqμ →∞ ) yields κ ℓ∞  =  148, 161, and 159 Wm−1K−1 for the 1NN, 2NN, and 3NN models, respectively. Thus, the extrapolation procedure makes the prediction more reliable and we achieve great agreements with experimental κℓ of the isotopically pure sample [38], 156 Wm−1K−1, especially for the 2NN and 3NN models. It should be noted, however, that the extrapolation procedure can overpredict experimental results because the phonon wave vector q and the mean-free-path ℓ(q, n)  =  |v (q, n)|τ (q, n) should be restricted in actual samples with grain boundaries and finite dimensions. In figure 4, we compare our computational results with experimental values and previous computational results. Here, we show our BTE-RTA results extrapolated to the infinite q density. In our calculations, phonon-isotope scatterings are not considered. Therefore, we compare the results with the experimental κℓ of isotopically enriched samples measured by Kremer et al [38]. Although an experimental κℓ above 300  K is not available for isotopically enriched systems, we can substitute the experimental result for the naturally occurring isotope concentrations, which is also shown in the figure, because phonon-phonon interactions dominate the

phonon scattering at high temperatures. As seen in the ­figure, our BTE-RTA result with the 2NN model agrees greatly with the experimental and previous computational results in a wide temperature range, verifying the accuracy of the extracted harmonic and anharmonic IFCs. In the high temperature range above 800 K, however, slight deviations from the experimental value can be seen. This cannot be explained by the effect of thermal expansion since a 0.5% increase in the lattice constant, which is slightly larger than the observed value at 1000 K [42], increased κℓ by about 4% . Therefore, the deviation implies the effects of the higher-order scattering process neglected here. In figure 4, we also show κℓ computed by NEMD simulations. In the following sections, we will present the details of our NEMD simulations and assess the quantitative reliability of the NEMD method. Since we will compare NEMD with our BTE results, we will solely employ the 1NN model for Si to reduce the computational costs. 3.3.  Non-equilibrium molecular dynamics

In NEMD simulations, the thermal conductivity is determined by Fourier's law: κ = −j / ∇T, (10)

where j is the heat flux vector and ∇ T is the temperature gradient. In figure 5, a simulation cell employed in this study is schematically shown. The system consists of two thermostat regions with different temperatures, central scattering region, and the fixed boundaries at each end. In the xy direction, the periodic boundary condition is used with the cross sectional area of 2.16 nm × 2.16 nm. In NEMD simulations, the lattice thermal conductivity can strongly depend on Lz, especially when the phonon mean-free-path (MFP) ℓ is larger than Lz. Therefore, in order to estimate bulk thermal conductivities, 7

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

3.4.  Size effects in NEMD

The errors of the bulk thermal conductivity in NEMD can be ascribed to the errors in κℓ(Lz) and those from the linear extrapolation (equation (11)). In general, κℓ(Lz) can be affected by the details of NEMD, such as the types and the length Lz of thermostats, the dimension in the xy direction, and the temperature difference ΔT. Among them, it has been reported that κℓ(Lz) considerably depends on LT, and small LT results in significant underestimation of the thermal conductivity [47]. Therefore, we employed much larger LT compared to the previous works so that long-wavelength phonons are efficiently equilibrated at each end. We also examined the dependence on the xy dimension by changing the cross sectional area from 2.16 nm × 2.16 nm to 3.24 nm × 3.24 nm, but did not observe significant differences. The size dependence in the z direction is two-fold and more complicated. Firstly, Lz limits the available phonon modes in the NEMD system, i.e. the phonon wavelength cannot be larger than Lz. Therefore, insufficient Lz causes an underestimation of κℓ as insufficient number of q points does in BTE. However, we consider such an effect is minor in NEMD simulations because Lz is far larger than the system size necessary for the convergence of κℓ in BTE (figure 3) and EMD-GK [11, 20]. Secondly, the system-size length Lz limits the MFP of phonons. If the MFP is independent on q and n, the linear relation (equation (11)) should be valid. However, Sellan et al [20] reported that the linear relation did not hold well in Si, because the phonon properties such as phonon MFP were considerably mode-dependent. They showed that the minimum Lz used for the extrapolation needed to be sufficiently large, or κℓ(∞) in NEMD could be underestimated compared to the results of BTE and EMD-GK. Therefore, we employed relatively large Lz up to 276 nm, even though that seemed to be insufficient for bulk Si. To elucidate the reason of underestimations, we compare the size-dependent thermal conductivity κℓ(Lz) in NEMD and that of BTE. Here, we employ the Matthiessen's rule to include the boundary effect in BTE (equation (9)), where the effective MFP of phonons is given by [48]

Figure 5. Schematic view of the supercell for NEMD simulations. The thermostat regions were attached to the central conduction region to generate heat flux.

κℓ(Lz) needs to be estimated for various Lz as explained below. In this study, the length of the scattering region Lz are 138, 207, and 276 nm (corresponding to 256, 384, and 512 conventional cell units). For generating non-equilibrium steady states, constant temperature methods and reversed cause-and-effect methods are available [43, 44]. In this study, we employed a constant temperature approach where the temperature in the heat source and the heat sink was controlled by Nosé–Hoover thermostats [45]. For thermostat regions, we employed LT  =  Lz/2, where individual Nosé variable was assigned to every 256 atoms. The heat flux was estimated by monitoring energy exchanges in the thermostat regions. To obtain smooth temperature profiles and the convergence of heat flux, NVT simulations were performed for 15–20 × 106 MD steps (corresponding to 15–20 ns). Upon integrating the equation of motion, the velocity Verlet method of Liouville-operator type was employed with the time step Δt = 1.0 fs, and reasonable conservation of the pseudo-Hamiltonian was confirmed. From the long-time simulation, we estimated the average and standard error of size-dependent thermal conductivity κℓ(Lz). For that purpose, we discarded first 5 ns of simulation data and uniformly divided the remaining data into segments. Each segment contains 1 ns simulation data, from which thermal conductivity of the segment is estimated. Then, the average and standard error were estimated following the usual statistical analysis as in reference [46]. To estimate the bulk thermal conductivity from NEMD simulations, we employ the linear extrapolation method proposed by Schelling et al [11] where the effective MFP of phonon ℓ eff is approximated as follows:

1 1 1 = + . (12) τeff ( q, n ) τ ( q, n ) τb ( q, n )

1 1 2 = + . (11) ℓ eff ℓ L z

We only consider the boundary along the z direction as in NEMD, where the boundary scattering probability is given by

Here, ℓ is the intrinsic phonon MFP and the second term is the effect of the boundary scattering at each side of the system. If phonon properties, such as the phonon velocity and the relaxation time, can be well represented by their average values, the thermal conductivity is proportional to the effective MFP and its size dependence is given by κ ℓ−1( L z ) = κ ℓ−1( ∞ ) + 2aL z−1 with a proportionality coefficient a. In figure 4, the lattice thermal conductivity of bulk Si obtained from the linear extrapolation is also shown. According to our estimate, NEMD considerably underestimated the lattice thermal conductivity of bulk Si. The magnitude of underestimation was larger in the low temperature region. In the following section, the reason of the underestimation is investigated by comparing with BTE.

α | v z ( q, n ) | 1 = . (13) τb ( q, n ) Lz

Here, we introduce the parameter α to explain the size dependence in NEMD. α  =  2 corresponds to the case where the MFP of phonons are restricted by the system size length Lz in NEMD simulations (see figure 5), which is consistent with equation (11). When α < 2, it means that the MFP is restricted by the characteristic length L′z = (2/α)Lz that is larger than Lz. Substituting τeff ( q, n ) for τ (q, n) in equation (9) yields ℓz ( q, n ) BTE ( L ) = 1 κ(14) . ∑ c ( q, n ) v z ( q, n ) z ℓ Ω Nq q, n 1 + α ℓz ( q, n ) / L z 8

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

thermoelectric material and its lattice thermal conductivity is much lower than that of Si. Since MFPs of heat-carrying acoustic phonons are expected to be small in Mg2Si, ℓ/Lz ratio should be readily reduced with moderate computational costs. Therefore, it is useful to investigate the nonlinear size-effect in NEMD and to assess the quantitative reliability of NEMD for systems with low κℓ. Firstly, we generated displacement-force data set by LDA calculations. For Mg2Si, we also employed a 3 × 3 × 3 conventional cell (324 atoms in total) to consider distant interactions. The cutoff energy of 33  Ry was employed with norm-conserving pseudopotentials. To determine harmonic IFCs, we employed 5  ×  5  ×  5 Monkhorst-Pack k grid to ensure the accuracy, whereas for anharmonic IFCs we performed FPMD simulations with 2 × 2 × 2 k grid. Since the temperature in FPMD is somewhat arbitrary, we employed 600  K in this study. Then, harmonic and anharmonic IFCs were extracted from the generated references. For harmonic terms, we considered Si–Si interactions up to fourth-NN, Si–Mg interactions up to fifth-NN, and Mg–Mg interactions up to seventh-NN shells. For cubic and quartic terms, we considered the nearest Si–Si, Mg–Mg, and Si–Mg interactions, whereas first-NN Si–Mg interactions were considered for fifth- and sixth-order IFCs. To reduce the computational costs, we only included two-body interactions for the quartic terms. The resulting number of IFCs were 70, 40, 60, 43, and 85 for second, third, fourth, fifth, and sixth orders, respectively. We found that the cubic IFCs extracted from the trajectory of FPMD agreed well with those extracted by the displacement method3. Finally, thermal conductivity calculations were performed both with BTE and NEMD. In BTE, only cubic anharmonicity was considered, whereas cubic and quartic terms were employed in NEMD. The details of computational conditions are basically the same as those of Si except cell dimensions. In NEMD simulations, cross sectional area of 2.51 nm × 2.51 nm was employed. To estimate bulk values, we employed Lz = 161, 241, and 321 nm (corresponding to 256, 384, and 512 conventional cell units). According to our analysis based on BTE, the effect of thermal expansion on κℓ was found to be as small as in bulk Si. Therefore, we performed BTE and NEMD simulations in a constant volume corresponding to the lattice constant at 0 K, 6.28 Å. In figure 7, we compare the computational results with experimental and previous computational results, which agree well with each other. Here, we show the BTE result extrapolated to the infinite q density. It should be noted that the isotope effect is not considered in this study, whereas the experimental data and previous computational results include such effect. Therefore, our prediction based on BTE gives slightly larger κℓ than the previous computational results. Li et al [51] reported that the phonon-isotope scattering effect reduced the room temperature κℓ of Mg2Si by 10%, which partly accounts for the difference. The experimental data of a polycrystalline sample shown in figure 7 includes contributions both from phonons and electrons. Due to the large band

Figure 6. Comparison of κ ℓ−1( L z ) versus L z−1 in BTE (solid and dashed lines) and NEMD (open symbols) simulations at 300 K. Here, the lattice thermal conductivity is normalized by κℓBTE obtained by equation (9) with 30 × 30 × 30 q grid. Linear extrapolation from the result of NEMD simulations is also shown (dotted lines).

Here, the MFP along the z direction ℓz(q, n) = |vz(q, n)|τ (q, n) is introduced, and the superscripts of κ are omitted since we focus on the heat transport along the z direction. In fi ­ gure 6, the calculated κ ℓ−1( L z ) at 300 K as a function of L z−1 is shown together with that obtained by NEMD simulations. To see the size dependence more clearly, we additionally performed NEMD simulations with smaller Lz. Here, we also show the result of the linear extrapolation indicated by dotted lines. As seen from the figure, κ ℓBTE ( L z ) with α  =  2 overestimate the boundary scattering probability of Si, whereas that with α = 1 well reproduces the size dependence of κℓ(Lz) in NEMD simulations at 300 K. For both α = 1 and α = 2, strong nonlinearities in κ ℓ−1( L z ) with respect to L z−1 are observed, in accord with the analysis by Sellan, showing the limitation of the linear extrapolation method (equation (11)). The linear extrapolation from three largest Lz in NEMD simulations yielded κℓ(∞) = 86 ± 4 Wm−1K−1, which is about 60% of our BTE result. If one can afford to employ larger Lz, a quantitative improvement is expected. Assuming that the size dependence in NEMD follows equation (14) and employing the linear extrapolation with Lz = L0 and Lz = 2L0, L0 ≳ 2700 nm is necessary to achieve 90% of the BTE result at 300 K. Therefore, it would be computationally expensive to achieve quantitative predictions of Si by NEMD. At high temperatures, much smaller L0 should be adequate because the MFP is decreased by stronger phonon-phonon scatterings. However, the discrepancy between BTE and NEMD is still noticeable at high temperatures. As increasing the temperature, the fourth- and higher-order anharmonicities absent in BTE become more important. Therefore, insufficient contributions from higherorder anharmonicities in BTE make κℓ larger. 3.5.  Applications to magnesium silicide

3

  (Supplementary data (stacks.iop.org/J.PhysCM/26/225402/mmedia)). Selected cubic IFCs of Mg2Si extracted by the displacement method and FPMD are provided online.

To see the validity of the computational approach, we also estimated the thermal conductivity of Mg2Si. Mg2Si is a 9

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

Figure 8. Mean-free-path of phonons in Si and Mg2Si at 300 K obtained by the perturbation method.

Figure 7. Temperature dependence of the lattice thermal conductivity of Mg2Si calculated by BTE-RTA (solid line) and NEMD (circle with error bars). The experimental values of total thermal conductivity [49] (triangle) and previous computational results by Chaput [50] (blue dotted line) and by Li et al [51] (green dashed line) are also shown.

modes, where each mode accounts for 26%, 35%, 35%, and 4% of the total κℓ at 300 K, respectively. If one can employ larger Lz (consequently smaller ℓ/Lz ratio), the severe underestimation can be partly cured as discussed in section 3.4. In figure 8, we show the MFP of phonons ℓ(q, n) = |v(q, n)|τ (q, n) for Si and Mg2Si at 300 K. As clearly seen in the figure, the MFP of phonons in Mg2Si is smaller than that in Si by one order of magnitude, which enables us to obtain the bulk thermal conductivity of Mg2Si by the linear extrapolation with Lz larger than 160 nm. In figure 6, we obtained good agreements between NEMD and κ ℓBTE ( L z ) with α  =  1. In the low temperature range, we observed that the κ ℓBTE ( L z ) with α  =  1 agreed better with NEMD than that with α  =  2. This implies the existence of a characteristic length L′z  =  2Lz that restrict the MFP of phonons. In our NEMD simulations, we took the thermostat regions large enough, i.e. LT = Lz/2, so that the phonon can be efficiently equilibrated at each end. Hence, the total length of the simulation cell L tot (excluding the fixed boundary) is L tot = L z + 2L T = 2L z, corresponding to L′z. In addition, when we employed thin thermostat layers, i.e. LT = 1.08 nm (containing 256 atoms) and L tot ≃ L z, we observed that the size dependence can be well explained by equation (14) with α ≃ 2. This indicates that the characteristic length L′z is manifested by L tot rather than Lz. To draw any conclusions about phonon scatterings at boundaries in NEMD, however, thorough study is necessary because it should depend on the detailed conditions of NEMD including kinds of thermostat. As increasing the temperature, agreements between NEMD and equation (14) become worse because of the higher-order anharmonicity in NEMD that is absent in BTE. For reliable and efficient computations of κℓ, it is advisable to select an appropriate method. In high-κℓ materials, such as bulk Si, fourth- and higher-order phonon scatterings are expected to be minor. Therefore, BTE should be suitable in this case because it can efficiently consider the long-wavelength phonons having long MFP that contribute significantly to heat transport. Applications of NEMD method to such high-κℓ systems, on the other hand, can lead to significant uncertainties in the lattice thermal conductivity due to

gap of Mg2Si, however, the electronic contribution should be negligible for the measured temperature range [49]. The agreement between NEMD and BTE is better in Mg2Si than in Si due to the small nonlinearity in κ ℓ−1( L z ) as shown in figure 6. The linear extrapolation from three NEMD data at 300  K yielded κℓ(∞)  =  13  ±  1  Wm−1K−1, which is larger than 90% of our BTE result. These results show the quantitative reliability of NEMD with first-principles force constants for Mg2Si. In figure 7, the results of NEMD seem to agree well with experimental and previous computational results. This accidental coincidence occurs due to the absence of isotope-phonon scatterings in our NEMD simulations. In this study, LO-TO splitting near the zone center is not considered for consistent comparison of BTE and NEMD. However, we assume that the error associated with neglecting the nonanalytic correction is minor because the heat is mainly carried by acoustic modes in simple materials as Mg2Si. 3.6. Discussions

The size dependence in κ ℓ−1( L z ) is quite different between Si and Mg2Si as shown in figure 6. This can happen because the size dependence should be affected by the detailed phonon properties. To analyze the strong nonlinearity in Si, we decomposed the thermal conductivity contributions into transverse acoustic (TA), longitudinal acoustic (LA), and optical modes. Then, we observed severe nonlinearities for acoustic modes, especially for the LA mode. Here, we evaluated the magnitude of nonlinearity for each mode by the ratio κ ℓ ( ∞ ) / κ ℓBTE. Since it is not possible to perform NEMD simulations for each mode, we alternatively estimated κ ℓ−1( L z ) by equation (14) with Lz = 138, 207, and 276 nm, from which κℓ(∞) was estimated by the linear extrapolation (equation (11)). The ratio of the extrapolated thermal conductivity to κ ℓBTE is 54% and 68% for two TA modes, 38% for the LA mode, and 98% for optical 10

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

NEMD can also be a powerful and reliable tool when applied to realistic systems with low thermal conductivities, because it can easily simulate the effects of impurities and interfaces, which should be of great significance in micro-/nanostructured devices consisting of thermoelectric materials.

the difficulty in removing the boundary effect by the linear extrapolation as discussed above. Fortunately, this difficulty is expected to be minor in low-κℓ systems such as thermoelectric materials and nanostructured devices with surfaces/ interfaces because the ℓ/Lz ratio can be readily reduced with moderate computational costs, as is the case for Mg2Si. An advantage of MD-based methods for such complex systems is more precise description of phonon scatterings with fourthand higher-order anharmonicities which are absent in BTE. In addition, compared to BTE, the computational efficiency of MD-based methods becomes better when applied to complex materials. The computational cost of BTE calculations is O ( Nc3Nq2 ) even if the lowest order perturbation is considered, where Nc is the number of atoms in the primitive cell. Therefore, it requires much computational power for complex thermoelectric materials with large Nc. The computational cost of MD simulations, on the other hand, scales linearly with system size because of the short-range nature of ALMs. Furthermore, MD methods can be applied to novel structures with double-well potentials as in inorganic clathrate compounds [52], to which BTE cannot be applied due to imaginary phonon branches. The applications to such a novel structure will be presented elsewhere.

Acknowledgments TT would like to thank A Togo for providing unpublished data of force constants in bulk Si. This study is supported by a Grant-in-Aid for Scientific Research on Innovative Areas ‘Materials Design through Computics: Complex Correlation and Non-Equilibrium Dynamics’ (grant no. 22104006) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. The computation in  this work has been done using the facilities of the Super­computer Center, Institute for Solid State Physics, The ­University of Tokyo. Appendix A.  Lifetime of phonons In calculating the relaxation time τ (q, n) in equation (9), we consider three-phonon scattering processes and approximate τ (q, n) as τ (q, n) ≈ 1/(2Γqn). Here, Γqn(ω) is the imaginary part of the phonon self-energy, which is given by [53]

4. Conclusions

Γqn ( ω ) =

We propose a method to extract anharmonic force constants from the trajectory of first-principles molecular dynamics simulations. Cubic, quartic, and higher-order force constants of bulk Si are calculated by numerical fitting and we obtain good agreements with the conventional finite-displacement method. Also, we consider the effect of temperature in firstprinciples MD on the anharmonic terms and find that the accuracy of cubic and quartic force constants are not sensitive to the temperature if higher-order corrections are included in the fitting. With this approach, cubic and quartic force constants, which characterize phonon-phonon scatterings in bulk system, can be readily calculated. The calculated harmonic and anharmonic force constants are combined with NEMD and BTE, and the lattice thermal conductivity of Si is estimated. The predicted values in BTE agree well with the experimental measurement except the high temperature region, which can be attributed to the absence of fourth- and higher-order phonon-phonon scatterings in BTE simulations. In NEMD, noticeable underestimations are observed in a wide temperature range. To understand the discrepancy, we directly compare the result of NEMD and BTE with boundary effects and observe the strong nonlinearity in size-dependent thermal conductivity, which is caused by acoustic phonons having long mean-free-path. Such a severe nonlinearity becomes noticeable in high-κℓ systems, which makes it difficult to estimate bulk values by the linear extrapolation method. We also calculate the lattice thermal conductivity of Mg2Si, a thermoelectric material. We observe that both BTE and NEMD agree well with experimental values and demonstrate the quantitative reliability of NEMD simulations for systems with low κℓ. On the basis of these results, we consider that

π 2

∑ ∑|

V−(3)qn, q1n1, q2n2 |2

q1, q2 n1, n2

× [( f1 + f2 + 1) δ ( ω − ω1 − ω 2 ) − ( f1 − f2 ) δ ( ω − ω1 + ω 2 )



+ ( f1 − f2 ) δ ( ω + ω1 − ω 2 )] .

(A.1)

Here, f is the Bose-Einstein distribution function and the subscripts qini are replaced by i for simplicity (e.g. fq1n1 → f1). The coefficient of three-phonon interactions V q(3) n, q′n′, q′′n′′ is defined by

Vq(3)n,q n ,q n

′ ′ ′′ ′′

3

=(

ℏ 2 ) 2Nq

×



× ×

1 ωqnω q′n′ω q′′n′′

∑

′ ′′ ei(q . m1+ q . m2+ q . m3)

m1, m2 , m3

∑ ∑

α1, α 2, α3 μ1, μ2 , μ3

Φ mμ11,αμ12,,mμ32α2, m3α3

e(μq1n)α1e(μq2′n′)α e(μq3′′n′′)α 2

Mα1Mα2Mα3

3

,

(A.2)

where Mα is the atomic mass and e(qn) is the polarization vector of phonons. To estimate equation (A.1), delta functions are usually 1 ϵ replaced by Lorentzian δ ( ω ) ≈ with smearing width π ω 2 + ϵ2 ϵ. Small ϵ is favorable to obtain detailed structures of twophonon density of states (DOS), but ϵ should be large enough to avoid numerical oscillations in the spectrum. Finding an appropriate value of ϵ is not a trivial task because it should depend on the q-grid density. In this report, we employ an 11

T Tadano et al

J. Phys.: Condens. Matter 26 (2014) 225402

[19] Esfarjani K and Stokes H T 2012 Phys. Rev. B 86 019904 [20] Sellan D P, Landry E S, Turney J E, McGaughey A J H and Amon C H 2010 Phys. Rev. B 81 214305 [21] Gonze X 1997 Phys. Rev. B 55 10337–54 [22] Baroni S, de Gironcoli S, Dal Corso A and Giannozzi P 2001 Rev. Mod. Phys. 73 515–62 [23] Parlinski K, Li Z Q and Kawazoe Y 1997 Phys. Rev. Lett. 78 4063–6 [24] Giannozzi P et al 2009 J. Phys.: Condens. Matter 21 395502 [25] Ceperley D M and Alder B J 1980 Phys. Rev. Lett. 45 566–9 [26] Perdew J P and Zunger A 1981 Phys. Rev. B 23 5048–79 [27] Monkhorst H J and Pack J D 1976 Phys. Rev. B 13 5188–92 [28] Berendsen H J C, Postma J P M, van Gunsteren W F, DiNola A and Haak J R 1984 J. Chem. Phys. 81 3684–90 [29] Nilsson G and Nelin G 1972 Phys. Rev. B 6 3777–86 [30] Kulda J, Strauch D, Pavone P and Ishii Y 1994 Phys. Rev. B 50 13347–54 [31] Aouissi M, Hamdi I, Meskini N and Qteish A 2006 Phys. Rev. B 74 054302 [32] Stillinger F H and Weber T A 1985 Phys. Rev. B 31 5262–71 [33] Dekura H, Tsuchiya T and Tsuchiya J 2013 Phys. Rev. Lett. 110 025904 [34] Turney J E, Landry E S, McGaughey A J H and Amon C H 2009 Phys. Rev. B 79 064301 [35] He Y, Savic I, Donadio D and Galli G 2012 Phys. Chem. Chem. Phys. 14 16209–22 [36] Tersoff J 1988 Phys. Rev. B 38 9902–5 [37] Glassbrenner C J and Slack G A 1964 Phys. Rev. 134 A1058–69 [38] Kremer R, Graf K, Cardona M, Devyatykh G, Gusev A, Gibin A, Inyushkin A, Taldenkov A and Pohl H J 2004 Solid State Commun. 131 499–503 [39] Esfarjani K, Chen G and Stokes H T 2011 Phys. Rev. B 84 085204 [40] Blöchl P E, Jepsen O and Andersen O K 1994 Phys. Rev. B 49 16223–33 [41] Lindsay L, Broido D A and Reinecke T L 2013 Phys. Rev. B 87 165201 [42] Okada Y and Tokumaru Y 1984 J. Appl. Phys. 56 314–20 [43] Ikeshoji T and Hafskjold B 1994 Mol. Phys. 81 251–61 [44] Müller-Plathe F 1997 J. Chem. Phys. 106 6082–5 [45] Nosé S 1984 J. Chem. Phys. 81 511–19 [46] Zhou X W, Aubry S, Jones R E, Greenstein A and Schelling P K 2009 Phys. Rev. B 79 115201 [47] Shiomi J and Maruyama S 2008 Japan. J. Appl. Phys. 47 2005–9 [48] Ziman J M 2001 Electrons and Phonons (New York: Oxford University Press) [49] Martin J J 1972 J. Phys. Chem. Solids 33 1139–48 [50] Chaput L 2013 Phys. Rev. Lett. 110 265506 [51] Li W, Lindsay L, Broido D A, Stewart D A and Mingo N 2012 Phys. Rev. B 86 174307 [52] Suekuni K, Takasu Y, Hasegawa T, Ogita N, Udagawa M, Avila M A and Takabatake T 2010 Phys. Rev. B 81 205207 [53] Maradudin A A and Fein A E 1962 Phys. Rev. 128 2589–608

alternative way to estimate equation (A.1). Due to the momentum conservation, equation (A.2) vanishes unless q+q′+q′′ is a integral multiple of the reciprocal lattice vector G. Therefore, the double sum over q points in equation (A.1) can be reduced to a single summation. For example, the first term of equation (A.1) can be written as π ∑ ∑ ∣ V−(3)qn, q1n1, q2′ n2 ∣2 × ( f1 + f2 ′ + 1) δ ( ω − ω1− ω2′ ) , (A.3) 2 n ,n q 1

2

1

with q′2 =  q−q1+G and f2 ′ = fq2 n2. Since q, q1, and q′2 are ′ restricted to be from the first Brillouin zone, q′2 can be uniquely determined. The summation over q1 in equation (A.3) is nothing but a weighted-DOS, hence can be estimated by the tetrahedron method [40]. References [1] Slack G A 1995 New materials and performance limits for thermoelectric cooling CRC Handbook of Thermoelectrics ed D M Rowe (Boca Raton, FL: CRC) pp 407–40 [2] Hicks L D and Dresselhaus M S 1993 Phys. Rev. B 47 16631–4 [3] Hicks L D and Dresselhaus M S 1993 Phys. Rev. B 47 12727–31 [4] Venkatasubramanian R, Siivola E, Colpitts T and O'Quinn B 2001 Nature 413 597–602 [5] Boukai A I, Bunimovich Y, Tahir-Kheli J, Yu J K, Goddard W A and Heath J R 2008 Nature 451 168–71 [6] Hochbaum A I, Chen R, Delgado R D, Liang W, Garnett E C, Najarian M, Majumdar A and Yang P 2008 Nature 451 163–7 [7] Broido D A, Ward A and Mingo N 2005 Phys. Rev. B 72 014308 [8] Dong J, Sankey O F and Myles C W 2001 Phys. Rev. Lett. 86 2361–4 [9] Donadio D and Galli G 2009 Phys. Rev. Lett. 102 195901 [10] He Y and Galli G 2012 Phys. Rev. Lett. 108 215901 [11] Schelling P K, Phillpot S R and Keblinski P 2002 Phys. Rev. B 65 144306 [12] Broido D A, Malorny M, Birner G, Mingo N and Stewart D A 2007 Appl. Phys. Lett. 91 231922 [13] Shiga T, Shiomi J, Ma J, Delaire O, Radzynski T, Lusakowski A, Esfarjani K and Chen G 2012 Phys. Rev. B 85 155203 [14] Shiomi J, Esfarjani K and Chen G 2011 Phys. Rev. B 84 104302 [15] Stackhouse S, Stixrude L and Karki B B 2010 Phys. Rev. Lett. 104 208501 [16] Huang B L and Kaviany M 2008 Phys. Rev. B 77 125209 [17] Lee Y and Hwang G S 2012 Phys. Rev. B 85 125204 [18] Esfarjani K and Stokes H T 2008 Phys. Rev. B 77 144112

12