Thermal conductivity modeling of hybrid organic-inorganic crystals and superlattices

Nano Energy 41 (2017) 394–407

Contents lists available at ScienceDirect

Nano Energy journal homepage: www.elsevier.com/locate/nanoen

Review

Thermal conductivity modeling of hybrid organic-inorganic crystals and superlattices Xin Qiana, Xiaokun Gub, Ronggui Yanga,c, a b c

MARK

⁎

Department of Mechanical Engineering, University of Colorado, Boulder, CO 80309, USA Institute of Engineering Thermophysics, Shanghai Jiao Tong University, Shanghai 200240, China Materials Science and Engineering Program, University of Colorado, Boulder, CO 80309, USA

A R T I C L E I N F O

A B S T R A C T

Keywords: Thermal conductivity Phonon transport Organic-inorganic Hybrid material Superlattice Organometal halide perovskite

Hybrid organic-inorganic materials have attracted intensive research interests due to the superb physical properties that take advantage of both organic and inorganic components and their promising applications in functional devices such as light emitting diodes, solar cells, and flexible thermoelectrics. Understanding thermal transport in hybrid materials is of importance for both the reliability and the performance of the hybrid materialbased systems. The organic-inorganic hybrid materials can be classified as nanocomposites, superlattices and crystals depending on the strength of organic-inorganic bonding and the feature size of organic and inorganic components. Previous research on thermal transport mainly focus on hybrid nanocomposites where phononinterface scattering dominates the thermal transport process. However, in hybrid superlattices and crystals, organic and inorganic components are blended at the atomic scale by chemical bonds without a clear interface, and the phonon transport physics is not well understood. This article reviews the recent progress in ab-initio modeling on phonon transport in hybrid crystals and superlattices. The predicted thermal conductivity of II–VI based organic-inorganic semiconductors, organometal perovskites and organic-intercalated TiS2 superlattices are reviewed and compared with experimental data when available. This review provides guidance for modeling thermal conductivity of materials with complex atomic structures from the first principles, which could be a good tutorial for both graduate students and experienced researchers who are interested in thermal energy transport in emerging hybrid materials.

1. Introduction Hybrid organic-inorganic materials (HOIMs) are synthesized by combining the organic and inorganic components over length scales from a few angstroms to ~ 100 nm [1,2]. The most appealing feature of HOIMs is that they can be conveniently designed for properties that are difficult to achieve in either organic or inorganic materials alone. For example, superb electronic properties can be derived from the inorganic component while still maintaining structural flexibility due to the blending with organic components [2]. In the past decade, a variety of novel HOIMs have been synthesized with promising applications like light emitting diodes [3,4], flexible solar cells [5–11], and flexible thermoelectrics [12,13], to name a few. However, thermal transport in HOIMs is not well understood although thermal conductivity greatly affects the thermal stability and the performance of the hybrid materialbased devices [14]. Due to the structural diversity and the wide range of organic-inorganic coupling strength in HOIMs, it would be helpful to categorize the HOIMs into different groups for a better understanding of

⁎

their mechanical and thermal properties. In Fig. 1, we illustrate that the organic-inorganic hybrid materials can be categorized into three groups (hybrid organic-inorganic nanocomposites, superlattices and crystals) according to the organic-inorganic coupling strength and the feature size of organic and inorganic components. The first group of HOIM is the hybrid organic-inorganic composites whose organic and inorganic components are coupled by weak van der Waals forces or hydrogen bonds. While organic-inorganic composite materials have been studied for many decades [15–18], significant efforts have been devoted to understanding the thermal and mechanical properties of nanocomposites when the feature size of the organic/inorganic structural units in these composites is on the order of a few nanometers. In these hybrid nanocomposites with distinguishable organic-inorganic interfaces, such as carbon nanotube-polymer composites [19,20], self-assembled materials [21–25] and nanocrystal arrays [26–28] and super-atomic crystals [29], thermal transport is determined by organic component, inorganic component, and most dominantly the phonon-interface scattering [27]. However, the organic-inorganic interfaces are not

Corresponding author at: Department of Mechanical Engineering, University of Colorado, Boulder, CO 80309, USA. E-mail address: [email protected] (R. Yang).

http://dx.doi.org/10.1016/j.nanoen.2017.09.047 Received 3 August 2017; Received in revised form 17 September 2017; Accepted 25 September 2017 Available online 28 September 2017 2211-2855/ © 2017 Elsevier Ltd. All rights reserved.

Nano Energy 41 (2017) 394–407

X. Qian et al.

Fig. 1. Schematic illustrating that hybrid organicinorganic materials can be classified into three groups according the bonding strength between organic and inorganic components and their feature sizes.

including organic-inorganic II–VI based hybrid crystals, the organometal halide perovskites, and the organic-intercalated TiS2 superlattice. This review provides a good tutorial for those who are interested in thermal energy transport in emerging hybrid materials.

present in hybrid organic-inorganic crystals, when the blending of organic and inorganic constituents happens at the atomic scale [30]. In contrast to the hybrid nanocomposites with either amorphous organic or inorganic structures, the hybrid organic-inorganic crystals possess long range periodicity. A few examples of hybrid crystals are the family of group II–VI element based hybrid crystals [31] and the hybrid organometal halide perovskites [11, 32–35]. The dominant coupling mechanisms between the organic and inorganic constituents in these hybrid crystals are either covalent bonds (II–VI based hybrid crystals) or ionic bonds (hybrid perovskites). The third group of HOIMs is the hybrid organic-inorganic superlattices, which lies between the hybrid crystals and the hybrid nanocomposites, in terms of both the length scale of feature size and the bonding strength between the organic and inorganic components, as shown in Fig. 1. The most notable examples of hybrid organic-inorganic superlattices are the atomic/molecular layer deposited (ALD/MLD) organic-inorganic superlattice [36,37] and the organic intercalated superlattice [12,13]. In these superlattices, the organic-inorganic coupling strength can be in a wide range, from the strong covalent bonding in ALD/MLD superlattices to the electrostatic forces in the organic-intercalated superlattices. The feature size of the organic and inorganic component can also vary from sub-nanometer to a few nanometers. Although a few experimental research found that the hybrid crystals and superlattices usually have low thermal conductivity [30,36], it is still not well understood how phonons transfer their energy in these hybrid crystals and superlattices, and how the organic components affect the thermal transport. In this Review, based on the authors’ extensive experiences on modeling the thermal conductivity of organic-inorganic II–VI based hybrid crystals, the organometal halide perovskites, and the organicintercalated TiS2 superlattice, deriving from the contemporary study of nanoscale heat transfer [38–43], a general first-principle-driven computational methodology for calculating the lattice thermal conductivity of complex materials is discussed. It is noted that the hybrid crystals and superlattices described in this Review are semi-conductive, and the dominant heat carriers are lattice vibration (phonons). In Section 2, we outline the computational challenges in modeling the thermal conductivity of hybrid crystals and superlattices. In Section 3, a simulation strategy integrating first principles calculation and molecular dynamics simulations is reviewed. Section 4 reviews the implementation and the application of this general strategy to the materials we recently studied,

2. Computational challenges Much progress has been made on modeling the thermal conductivity of nanostructured materials over the past two decades [44–49]. The most distinguished examples are first-principles-based simulation methods like Boltzmann transport equation (BTE) or ab initio molecular dynamics (AIMD), taking great advantage of the progress in computational power [50–53]. Without fitting parameters, such first-principlesbased methods has been used to model the thermal conductivity of inorganic bulk crystals such as diamond [54,55], silicon [50,52,56], GaN [51] and two-dimensional materials such as graphene [38,57–59], silicene [38,60], black phosphorus [61–64] and many transition metal dichalcogenides [39,65,66]. Not only great agreements between theoretical calculations and measurement results on the thermal conductivity of materials have been found, but also these tools serve significantly speed up the materials discovery in both extremely high/low thermal conductivity and multifunctional thermal materials [50, 52–56, 66–67]. However, much less effort has been devoted to study phonon transport in hybrid organic-inorganic crystals and superlattices, not catching up the pace on the synthesis efforts for multifunctional thermal materials. This is understandable because there are significant computational challenges due to the structural complexity of hybrid organic-inorganic crystals and superlattices. For example, hybrid crystals usually have very large unit cells containing dozens of atoms, in comparison with only one to a few atoms in inorganic crystals. In addition the organic components usually have internal degrees of freedom of motion (for example, the organic ions in hybrid perovskite can rotate) instead of just vibrating around the local equilibrium [68–70]. To capture such dynamical disorder caused by the internal degrees of freedom of motion, a very large supercell is required to simulate phonon transport instead of just one unit cell with periodic boundary conditions. In superlattices, the organic components might even be amorphous. The lack of periodicity in the hybrid superlattices and the dynamical disorders in the hybrid crystals hinder the implementation of 395

Nano Energy 41 (2017) 394–407

X. Qian et al.

first-principles-based simulation methods like BTE or AIMD. In the firstprinciples-based BTE approach, thermal conductivity is calculated using phonon properties including phonon dispersion and detailed phonon-phonon scattering rates. However, the computation time increases dramatically with the number of atoms in a unit cell, making it a great challenge for hybrid materials. Similarly, the computational time required for direct AIMD simulations is not affordable to capture the size effect on the thermal conductivity of hybrid materials [71]. Alternatively, classical molecular dynamics (MD) simulation which traces atom movements according to Newton's second law of motion are computationally efficient [72]. However, it suffers from the lack of appropriate potential fields that govern interatomic interactions. Most of the potential fields in literature are not designed for hybrid organicinorganic crystals and superlattices. For example, the most commonly used potential fields including the Tersoff potential [73–75], StilingerWeber potential [73,76] and bond order potential [77–80] were developed for inorganic materials. However, the structure of the inorganic component in the hybrid materials can be very different from the original inorganic structure as some of the bonds are now connected to organic components, making it inappropriate to directly apply the potential fields in literature [81]. Another challenge is the interaction between the organic and the inorganic components, although the organic components can be reasonably well described using the potential fields such as the class II force field [82,83] and consistent valence force field (CVFF) [84]. To address these challenges, we present in this Review a strategy for calculating thermal conductivity of hybrid crystals

and superlattices using the ab-initio based simulations. The predicted thermal conductivity of II-VI based organic-inorganic semiconductors, organometal perovskites and organic-intercalated TiS2 superlattices are reviewed and compared with experimental data when available. 3. Strategy for ab-initio based simulation of thermal conductivity Fig. 2 summarizes the workflow for modeling the thermal conductivity of organic-inorganic hybrid crystals and superlattices by integrating the first-principles density functional theory (DFT) calculations and the classical molecular dynamics (MD) simulations. First-principle calculations are used to construct the potential fields. The very first step to construct the potential field is to assign functional forms. In HOIMs, the bonding nature of the organic and inorganic parts could be very different. Therefore, the potential field EV is separated into three parts: the potential for the organic part E O , the potential for the inorganic part EI and the organic-inorganic coupling EI − O :

EV = EI + E O + EI − O

(1)

Based on the bonding characteristics, a functional form is then assigned. For example, the Lennard-Jones potential is assigned to a van der Waals bond. Such a separation of the potential into three parts helps to greatly reduce the amount of work for developing the potential fields, by allowing to use some of the existing potential fields as parts for EI , E o , and EI − O . For example, the class II force fields [82,83] or the consistent valance force field (CVFF) [84] can be used to describe the Fig. 2. Computation strategy for obtaining thermal conductivity and phonon properties using first-principles based atomistic simulations.

396

Nano Energy 41 (2017) 394–407

X. Qian et al.

covalent bonding within the organic ligands. After assigning the functional forms to a potential field, a set of unknown parameters P in the functional form need to be determined. The parameters P essentially control the shape of the empirical potential surface and they are obtained by minimizing the difference between the empirical potential surface EV and the DFT energy surface EDFT :

W (P ) =

∑

EDFT ({ri}) − EV (P , {ri})

kα =

(2)

using optimization algorithms like quasi-Newton algorithm [85] or genetic algorithm [86], where ri denotes the atomic coordinate configurations. The unknown parameters P can also be obtained by the force matching method. Instead of the static sampling of energy surface, an ab-initio molecular dynamics (AIMD) simulation is performed, and the atomic trajectory ri (t ) and the interatomic forces FAIMD (ri (t )) are recorded. For any set of parameters, the empirical potential field would also predict a set of interatomic forces FV (P; ri ) = −∇ri EV (P; ri ) . Then the unknown parameters P can be obtained by minimizing the difference between the atomic forces from AIMD calculation and that calculated from the empirical potential field:

W (P ) =

1 Nt N

Nt

N

∑∑ t=1 i=1

FV (P; ri (t )) − FAIMD (ri (t )) 2

∫0

∞

Jα (0)⋅Jα (t ) dt

(4)

where α denotes the directions in the Cartesian coordinates, Jα (0)⋅Jα (t ) is the heat current autocorrelation function (HCACF), V is the volume of the simulation cell, kB is the Boltzmann constant and T is the temperature. A great advantage of EMD is that the thermal conductivity along different directions can be calculated in a single EMD run, and therefore EMD is suggested for calculating the thermal conductivity of hybrid crystals and superlattices. Noting that extracting thermal conductivity from MD simulations requires a series of convergence tests with respect to simulation cell size, integration time, etc. Readers who are interested in MD simulations for thermal conductivity are referred to literature [40,41,88,91,94–96]. In addition to the thermal conductivity, the MD simulation can also provide phonon properties, by projecting the atomic movements in real space into the momentum-frequency space where the phonons transport and scatter with each other. The most commonly used phonon modal analysis based on MD simulations include the normal mode analysis (NMA) [49,96] and spectral energy density (SED) [97,98]. The NMA projects the atomic vibrations into normal coordinates, and the phonon lifetime is extracted from the autocorrelation of the modal vibrational energy. Therefore, NMA requires eigenvectors as input which are usually determined from lattice dynamics calculation at 0 K. However, in some materials with strong anharmonicity, the eigenvector (phonon dispersion) at high temperature would shift away significantly from the 0 K lattice dynamics calculation [97]. The advantage of SED is that this method does not require predetermined eigenvectors. The SED analysis directly calculates the vibrational energy distribution in the frequency-momentum space. The SED ϕ (κ , ω) is defined as:

2

{ri}

V 3kB T 2

(3)

In most cases, both the energy surface fitting by Eq. (3) and the force matching using Eq. (4) can be used to obtain the parameters with subtle difference, but the force matching method is suggested when there are internal degrees of freedom. One example is the hybrid perovskite CH3NH3PbI3 (MAPbI3) whose organic ions CH3NH3 (MA) can rotate at high temperature, therefore a local sampling by small displacements would fail at large rotational degrees and overestimate the rotational energy barrier (see detailed discussion in next section). After parametrizing the potential field, a validation process need to be performed before calculating the thermal conductivity. To ensure proper fitting of the potential fields is achieved, the physical properties such as lattice constants and elastic properties that can be more easily calculated from both MD simulations and DFT calculations are often calculated and compared, or even compared with measurements when available. After a satisfying potential field is found, the thermal conductivity can be calculated rather straightforwardly using the classical molecular dynamics. The most commonly used classical MD methods to derive thermal conductivity can be classified into the non-equilibrium molecular dynamics (NEMD) [42,87–89] and equilibrium molecular dynamics (EMD) [41,90–91]. One of the mostly used NEMD method is the direct method, which uses a hot reservoir and a cold reservoir to create a steady state temperature gradient ∇T , and the heat flux q is calculated. Then the thermal conductivity is calculated by Fourier's law: k = − q/∇T , where ⋅ denotes the ensemble average. The direct method, however, usually requires a large simulation cell to establish a reproducible diffusive temperature profile [88]. In addition, the direct method only calculates thermal conductivity along the temperature gradient, making the computation costly when the thermal conductivity of a material is anisotropic, as multiple simulation runs need to be performed for thermal conductivity in different directions. Similar to the laser flash analysis measurement [92], a transient method called the approach-to-equilibrium molecular dynamics (AEMD) [89] uses a pulsed heating reservoir and monitors the temperature decay as a function of time. The thermal conductivity is then extracted using the characteristic decay time. This method, however, is also limited to isotropic materials. Different from the direct method and the AEMD method which drive the system out of thermal equilibrium, EMD calculates the thermal conductivity at thermal equilibrium using the Green – Kubo relation [93]:

2

ϕ (κ , ω) =

1 4πτ0 NT

τ0

∑ ∑ mb ∫0 ∑ vα ⎛⎜bl , t ⎞⎟⋅exp ⎡⎢iκ⋅r ⎛⎜ 0l ⎞⎟ − iωt⎤⎥ dt α

b

l

⎝

⎠

⎣

⎝ ⎠

⎦

(5) where κ is wavevector and ω is the phonon frequency, NT is the number

l of unit cells, vα ⎛⎜ ,t ⎞⎟ denotes the velocity component of the b -th basis ⎝b ⎠ atom in the l -th unit cell along the α -th direction, r ⎛⎜ l ⎞⎟ denotes the 0 ⎝ ⎠ position of the l -th unit cell, mb is the atomic mass of the b -th basis atom in the unit cell. In the (κ , ω) space, the location of SED peaks (bright colors) indicates the phonon dispersion and the bandwidth of the branches are related to the phonon lifetime. The smaller the broadening of the SED peak, the longer the phonon lifetime [12,31]. The phonon lifetime can be extracted by fitting the SED peaks using the Lorentzian function ϕ (κ, ω) =

I 1+(ω − ω0)2 /Γ 2

(6)

where I is the peak magnitude, ω0 is the frequency at the peak center and Γ is the bandwidth. The phonon lifetime can therefore be obtained by τ (κ ,ω) = 1/2Γ [49]. In the next section, we discuss the implementation of this integrated DFT-MD simulation strategy for predicting the thermal conductivity of II-VI based hybrid crystals [81], the organometal halide perovskites [99] and the organic-intercalated TiS2 superlattice [12]. 4. Thermal conductivities of hybrid organic – inorganic crystals and superlattices In this section, we discuss the application of the integrated DFT – MD simulation strategy to model the thermal conductivity of II–VI based hybrid organic – inorganic crystals, hybrid organometal halide perovskites, and organic – intercalated TiS2 superlattices. 397

Nano Energy 41 (2017) 394–407

X. Qian et al.

Fig. 3. Atomic structures of the II–VI based hybrid organic-inorganic crystals. The three black arrows denote the crystal orientations. The x-axis is along the direction where the organic-inorganic layers are stacked together. Parallel to the ZnTe monolayers, the y-axis is along the atomic ridges of the ZnTe monolayers, and z-axis is in the perpendicular direction. (a) The unit cell of β-ZnTe(en)0.5. (b) The unit cell of α-ZnTe(en)0.5. (c) From up to down: the folding of Zn-Te hexagonal rings along the short diagonals in the β-phase, the upper view of β -ZnTe monolayers and the side view of the β -ZnTe monolayers along the blue arrow. (d) From up to down: the folding of Zn-Te hexagonal rings along the longest diagonals in the α-phase, the upper view of α-ZnTe monolayers and the side view of the α-ZnTe monolayers along the blue arrow. (e) Atomic structure of ZnTe(pda)0.5. (f) Atomic structure of ZnTe (ptda)0.5.

organic component on the thermal conductivity of the hybrid crystals. As discussed in Section 3, the interatomic potential EV of α -ZnTe (en)0.5 can be divided into the three parts: the inorganic potential EI within the ZnTe monolayers, the organic potential E O within the organic ligand layers, and the coupling EI − O between the organic and inorganic constituents. Functional forms need to be assigned to different interaction terms. Appropriate choice of functional forms is essential for the success of this strategy. For thermal conductivity calculations of solids with strong chemical bonds, the potential does not need to be reactive (i.e. capable of describing the bond breaking process) because the atoms usually vibrates locally around the equilibrium position [100,101]. Therefore, simple functional forms with minimal number of empirical parameters are preferred for faster computation. For α-ZnTe(en)0.5, the potential EI , E O and EI − O can be further divided into different terms accounting for pairwise, angular and dihedral deformations. A list of typical deformations included in the potential field are shown in Fig. 4. For example, the inorganic potential energy can be divided into pairwise interaction EPI and angular interaction EAI :

4.1. II–VI based hybrid organic-inorganic crystals The II–VI based hybrid organic-inorganic crystal family has the chemical formula XY(L)0.5 where X is a group II-B element (Zn, Cd and Hg), Y is a group VI-A element (S, Se and Te), and L denotes the organic ligand. The II–VI based hybrid crystals form alternating organic-inorganic layered structures, which is a quantum-well designed for white light emitting diodes applications [3]. Among this family, the β-ZnTe (en)0.5 is the first discovered member in 2000 as the first example of crystalline hybrid organic-inorganic network formed by strong covalent bonds [31]. The hybrid crystal β-ZnTe(en)0.5 has a layered structure of alternating ZnTe monolayers connected by the organic ethylenediamine (CH2NH2)2, hereafter denoted as (en) ligands through covalent bonds as shown Fig. 3a. Besides the β-ZnTe(en)0.5, there are also other polymorphs including the α-ZnTe(en)0.5 (Fig. 3b) whose inorganic ZnTe monolayers has a different structure from the β- phase. Different from the β- phase in which the ZnTe hexagons are folded along the shorter diagonals (Fig. 3c), the α-ZnTe monolayers are formed by folding the ZnTe hexagons along the longest diagonals (Fig. 3d). Polymorphs like ZnTe(pda)0.5 and ZnTe(ptda)0.5 can form by replacing the (en) ligand in α-ZnTe(en)0.5 to longer organic molecules like H2N-(CH2)3-NH2 (pda) and H2N-(CH2)5-NH2 (ptda), as shown in Fig. 3e–f. Details of thermal conductivity modeling in the β-ZnTe(en)0.5 can be found in Ref. [81]. In this part, α-ZnTe(en)0.5 is first discussed in detail as an example to review the implementation of the integrated DFT-MD simulation methodology, then the thermal conductivity of other polymorphs with longer organic ligands like ZnTe(pda)0.5 and ZnTe(ptda)0.5 which were not published before, are discussed here to understand the effect of

EI = EPI + EAI =

∑ ϵ Ip (rp) + ∑ ϵaI (θa) a

p I

EP

I

EA

(7)

where ϵ Ip (rp) is the bonding energy of the atomic pair p (e.g. Zn-Te), depending on the atomic distance of the atomic pair rp , and ϵ aI (θa) is energy of an angle θa formed by an atomic triplet (e.g. Zn-Te-Zn). Morse potential is chosen to describe the pairwise energy of nearest neighbors (Zn–Te bonds) and the second nearest neighbors within the inorganic part: 398

Nano Energy 41 (2017) 394–407

X. Qian et al.

Fig. 4. Construction of the potential field for α-ZnTe(en)0.5. The potential field is divided into three parts: inorganic potential EI , organic potential E O and the organic-inorganic coupling EI − O . For each part, the potential can be further divided into different interaction terms to describe bond stretching, angular bending, dihedral bending, and the van der Waals coupling.

ϵ Ip (r ) = D (exp[−αm (r − r0)]−1)2

In Eq. (11), EPI− O includes the Zn-N bonding and the second nearest neighbor interactions (Zn-C and Te-N) and EAI− O includes the angular deformation energies of the triplets Zn-N-C, Zn-N-H and Te-Zn-N. Similarly, Morse potential and cosine-squared potential are used for EPI− O and EAI− O , respectively. After the assigning the functional forms to the potential, a set of unknown empirical parameters P in the functional form is then obtained by fitting the energy surface obtained from the DFT calculation. Fig. 5 shows an example of the DFT energy surface sampling and fitting process. The DFT energy surface sampled by displacing a Te atom along the Zn–Te bond. By varying the parameters P , the shape of the empirical potential surface would also change accordingly. Quasi-Newton algorithm [85] is used to adjust the parameters P until the difference between the DFT energy surface and the empirical potential surface reaches minimum. The parameters P at the best fit is then adopted for the optimized potential field. The developed potential field is then validated by comparing the elastic constants that can be easily calculated by both DFT simulations and MD simulations using the empirical potential field. The elastic constants are second-derivatives of the potential surface against dif1 ∂2E ferent deformations Cij = V ∂e ∂e , where ei and ej are components of the

(8)

where r0 is the bonding length at equilibrium. The characteristic bonding energy D and the characteristic potential well width αm are the empirical parameters to be determined. For three-body angular interaction, a cosine-squared form is used:

ϵ aI (θ) = K (cos θ − cos θ0)2

(9)

where the equilibrium angle θ0 and the angular stiffness K are the empirical parameter. Morse potential and the cosine-squared angular potential has been proven accurate and is commonly used to describe other covalent bonded systems like metal-organic-framework [102], bismuth telluride [103] and skutterudites [104]. For the organic part E O , the COMPASS potential [83] is directly applied for the pairwise, angular and dihedral deformations within each organic molecule. The inter-molecular coupling is described by the Lennard-Jones potential [105]:

σ 12 σ 6 ϵLJ = 4ϵ ⎛ ⎛ ⎞ − ⎛ ⎞ ⎞ ⎝r ⎠ ⎠ ⎝⎝ r ⎠ ⎜

⎟

(10)

i

Finally, the organic-inorganic coupling potential is also separated into pairwise deformation energy EPI− O and angular deformation energy EPI− O :

EI − O = EPI− O + EAI− O

j

strain tensor in Voigt notation [106], and V is the volume of the simulation cell. Table 1 shows great agreement of the lattice constants and Fig. 6 shows the agreement of elastic constants between the calculations from MD using the empirical potential field and the DFT calculations.

EI − O

(11) 399

Nano Energy 41 (2017) 394–407

X. Qian et al.

Fig. 7. Comparison of elastic constants in α-ZnTe(en)0.5, ZnTe(pda)0.5 and ZnTe(ptda)0.5.

Fig. 5. DFT potential energy surface sampled by displacing the Te atom along the Zn–Te bond is fitted by the empirical potential field. The stretching of Zn–Te bond shown in the figure is not to scale.

major deformation occurs in the dihedrals including the C-C-C-C or C-CC-N, yet these dihedrals have very small stiffness. As a result, the crystals are more flexible to shearing deformations when the organic ligand becomes longer. Fig. 8 shows the anisotropic thermal conductivity of α-ZnTe(en)0.5, ZnTe(pda)0.5 and ZnTe(ptda)0.5. The length of the organic ligands affects thermal conductivity anisotropically. The thermal conductivity parallel to the ZnTe monolayers (k y and k z ) decreased when the organic chains becomes longer. Due to the decreased shearing stiffness C55 and C66 , the group velocity of TA phonons would also decrease along the yand z- axis, resulting in the reduced k y and k z (see Table 2). The thermal conductivity along the stacking direction k x , however, has no significant changes, because the increased LA group velocity and the decreased TA group velocity affect k x counteractively.

Table 1. Lattice parameters (Å) of α-ZnTe(en)0.5. Parameters Lattice constants

a b

DFT x y z

7.241 7.025 17.815

MD 7.230a 7.023a 17.715a

7.075 7.055 17.813

Experiment 7.061b 6.927b 17.524b

Ref. [107]. Ref. [108].

Using the similar method, the potential fields for polymorphs with longer organic chains ZnTe(pda)0.5 and ZnTe(ptda)0.5 are also developed. Fig. 7 compares the elastic constants of α-ZnTe(en)0.5, ZnTe (pda)0.5 and ZnTe(ptda)0.5. Interestingly, the normal stiffness C11 along the stacking direction (x-axis) increased when the length of the organic ligand increases. When the same normal strain is applied along the xaxis, more C-C-C triplets are involved in the deformation when the organic ligand becomes longer. Therefore, the energy change associated with the applied strain increases with the organic ligand length, resulting in a larger C11 for the polymorphs with longer organic ligands. However, the shearing stiffness C55 and C66 decreased dramatically due to the increased degrees of freedom of dihedral deformations in the organic ligands. When the crystal is imposed with shearing deformations parallel to the ZnTe monolayer (corresponding to C55 and C66 ), the

4.2. Organometal halide perovskites The organometal halide perovskite is a family of hybrid organicinorganic crystals with chemical formula XYZ3, where the X site is occupied by organic ions (e.g. CH3NH3, hereafter MA), the Y site is a metal element (e.g. Sn, Pb etc.) and the Z site is a halide element (Cl, Br and I) [110]. These hybrid perovskites have attracted intensive research interests in photovoltaics [9,11,111–114] and recently in thermoelectrics [115,116]. However, thermal stability remains a problem for photovoltaics based on these hybrid perovskites [117]. A phase change from the tetragonal phase to the pseudocubic phase would occur under

Fig. 6. (a) Elastic constants of α-ZnTe(en)0.5 calculated from MD and DFT simulations. (b) The strains corresponding to the elastic constants are indicate by arrows.

400

Nano Energy 41 (2017) 394–407

X. Qian et al.

potential well is sampled, the empirical potential surface would dramatically overestimate the rotational energy barrier (~ 0.54 eV by DFT). At 400 K, the thermal excitation energy ~ 0.414 eV (kB T times 12 atoms in a unit cell) is comparable with the rotational energy barrier, making the rotations of MA ions non-negligible. By matching the AIMD forces at 400 K, the rotational barrier can be reasonably represented by the empirical potential field (~ 0.57 eV). After the parameters P in the functional form are determined, the potential is verified by comparing the vibrational density of states from MD to the DFT calculation [120], see details in Ref. [102]. With the developed potential field, the thermal conductivity of MAPbI3 and MASnI3 are then calculated using the classical MD simulations, as shown in Fig. 11a. Despite the detailed curvature discrepancy of the rotational potential surface between the fitting and the DFT calculations, the thermal conductivity of MAPbI3 still shows good agreement with experiment. This is because the rotational modes of MA ions have negligible contributions to thermal conductivity [99]. However, accurately representing the rotational barrier is important. If the rotational barrier of MA ions is severely overestimated, then these rotational modes are constrained at room temperature. As a result, the phonon-phonon scattering rate involving the MA rotations would be underestimated and thereby the thermal conductivity would be overestimated [121]. Both the MAPbI3 and MASnI3 have very low thermal conductivity in the tetragonal phase, but the thermal conductivity of MAPbI3 (0.59 W/ m K) is slightly lower than MASnI3 (0.69 W/m K) at around 300 K. This is because the Sn–I bond is shorter and stiffer than the Pb–I bond, resulting in a higher acoustic velocity in MASnI3 [33]. In both materials, a sharp increase in the thermal conductivity is observed during the transition from the tetragonal phase to the pseudocubic phase. However, the thermal conductivity of the pseudocubic MASnI3 is smaller than MAPbI3. This is because the third order force constant in the Sn–I bond is ~ 20% larger than that in the Pb–I bond, which results in the stronger anharmonicity. With the lower thermal conductivity and the much smaller electric resistivity in the pseudocubic MASnI3 [116], it is expected that replacing the Pb atoms with Sn atoms would significantly improve the thermoelectric performance. Fig. 11b compares the computational results of thermal conductivity in MAPbI3 existed in the literature and the measurement by Pisoni el. al [117]. Among the computational studies, Hata et al. [121] used a similar ab-initio based force matching method, but their potential overestimate the results above 50 K. The possible reason is that all bonds are assumed harmonic in their simulations, which significantly underestimate the anharmonicity. Another noteworthy potential field for thermal conductivity modeling of MAPbI3 is the MYP potential with purely pairwise interaction [68], which is adopted for MD simulation by Wang et al. [122] and Caddeo et al. [123] This MYP potential field, however, has relatively low accuracy of thermal conductivity depending on the simulation method (EMD [122] or AEMD [123]). Considering MYP potential is primarily designed to reproduce atomic structures for each phase of MAPbI3 with one set of parameters, the thermal conductivity accuracy for each phase might be compromised. Recently the direct non-equilibrium ab-initio molecular dynamics (NEAIMD) simulation without any fitting parameters has been proposed by Yue et al. [124] to model thermal transport in the hybrid

Fig. 8. Comparison of thermal conductivity in α-ZnTe(en)0.5, ZnTe(pda)0.5 and ZnTe (ptda)0.5.

continuous solar radiation (see Fig. 9) [110]. A good understanding in thermal conductivity of hybrid perovskites is therefore of great importance. Different from the covalently bonded hybrid crystals as described in Section 4.1, the organic-inorganic coupling in MAPbI3 and MASnI3 is ionic, and the Buckingham potential and Coulombic potential are therefore used to describe EI − O [118]:

EI − O =

rij

∑ A exp ⎛− ρ ⎞ − i, j

⎜

⎟

⎝

⎠

qi qj C + 6 rij rij

(12)

where i is the index of an inorganic atom (Pb or I) and j is the index of an organic atom (C, N or H). A , ρ and C are unknown parameters to be determined, and q is the atomic charge. To accurately determine the Coulombic energy, Bader analysis based on DFT simulation is performed to calculate the atomic charges [119], see details in Ref. [99]. Functional forms for EI and E O are also given in Ref. [99]. However, the method of matching the local energy surface as described in part 4.1 becomes inappropriate for modeling the material systems with internal degrees of freedom, because the non-zero temperature would play an important role in affecting the rotational dynamics of MA ions. At low temperature, the atom can only “feel” the potential surface very close to the static equilibrium position, but at the elevated temperature the internal rotational degrees of freedom are introduced and the potential surface far away from the equilibrium position becomes important. The DFT calculations by imposing displacements to the MA ions which has large internal rotation become very large to ensure that the potential surface along different rotation paths are effectively sampled. To solve this problem, ab-initio molecular dynamics (AIMD) is used to sample the potential surface randomly by the random atomic vibration. The unknown parameters P of the potential field is then obtained by matching the AIMD forces for minimal W (P ) in Eq. (3). Fig. 10 shows the potential surface associated with the rotation of the MA ion around the c-axis. If only a small portion of the

Table 2. The group velocities (m/s) of ZnTe(en)0.5, ZnTe(pda)0.5 and ZnTe(ptda)0.5 by solving the anisotropic wave propagation equation using the Christoffel package [109] .

x

LA TA1 TA2

z

y

en

pda

ptda

en

pda

ptda

en

pda

ptda

3449.7 1242.5 1141.9

4537.3 1192.2 909.9

5406.5 622.1 568.8

2962.2 1547.6 1252.5

3322.8 1518.3 1192.2

3685.7 1737.4 622.1

3857.9 1547.6 1141.9

4055.7 1518.3 909.9

3941 1737.4 568.8

401

Nano Energy 41 (2017) 394–407

X. Qian et al.

Fig. 9. The crystal structure of (a) the tetragonal and (b) the pseudocubic phases of MAPbI3 (MASnI3) projected in (010) and (001) plane from left to right correspondingly. The a and b axes of the both phases are defined along the two directions with smaller lattice constants, and the c axis is defined along the direction with the largest lattice constant. The legend on the top indicates different elements depicted in the figure.

[99], the NEAIMD also predicted the increased thermal conductivity during the tetragonal-to-pseudocubic phase transition occurs in MAPbI3 at 330 K. By analyzing the phonon properties from AIMD, they observed that the modes above 25 THz are purely organic vibrations (Fig. 12a) with negligible group velocities. The phonon modes above 25 THz make small contribution to the thermal conductivity despite their relatively long lifetime. Phonons modes below 5 THz contribute to majority of the thermal conductivity, which are vibrations of Pb-I and movements of the MA ions as a whole. 4.3. Organic-intercalated TiS2 superlattices The successful synthesis of organic-intercalated TiS2 superlattice realized for the first time flexible n-type thermoelectric material with high thermoelectric figure of merit (ZT) [12,13]. The TiS2 intercalated with the organic hexylammonium ions were found to have only 1/30 the in-plane thermal conductivity (0.12 W/m K) of the pristine TiS2 (4.2 W/m K), leading to a high ZT = 0.28 at 373 K. To understand the origin of such a dramatic thermal conductivity reduction, the integrated DFT-MD method has been applied to analyze the phonon transport in organic-intercalated TiS2 superlattice. There is a significant difference between the organic-intercalated superlattices and the hybrid organic-inorganic crystals. In the hybrid crystals, the coupling between the organic and inorganic parts is so strong that the hybrid potential field cannot be determined separately from the isolated organic and inorganic components. For example, the

Fig. 10. Comparison between the energy surfaces from DFT calculations and the empirical potential which fits the local curvature of the DFT energy surface. The rotational axis of MA ion is in the [001] direction.

organic-inorganic perovskites. However, NEAIMD can only be used to model a system with very small number of atoms (a few hundred) due to the high computational cost, which limits its application to materials with low thermal conductivity. Similar to Qian's MD simulation results

402

Nano Energy 41 (2017) 394–407

X. Qian et al.

inorganic Pb-I cages in the perovskite MAPbI3 are not stable without the organic ions. Therefore, the potential EI cannot be determined without considering the effect of organic ions. The situation is different in the organic-inorganic superlattices. In the organic-intercalated TiS2 superlattice, the organic intercalants only expand the van der Waals gaps without breaking the chemical bonds within the TiS2 monolayers. Therefore, the potential for the inorganic TiS2 can be modeled separately based on bulk TiS2. In layered TiS2, EI can be separated into two parts: I I EI = Ecov + ELJ

(13)

I Ecov

where is the covalent bonding energy within each TiS2 layer (taking I is the the advantage of potential for the pristine TiS2 crystal), and ELJ inter-layer coupling described by Lennard-Jones potential. Taylor expansion up to the fourth order is used to approximate the covalent bonding in TiS2: I Ecov = E0 +

1 + 4!

1 2!

∑ ∑ ϕijαβ ujβ + ij

∑∑

αβ

1 3!

∑ ∑ ψijkαβγ uiα ujβ ukγ ijk αβγ

αβγδ α β γ δ χijkl ui uj uk ul

(14)

ijkl αβγδ

where E0 is the energy when all atoms are sitting in their equilibrium positions, uiα , ( ujβ , ukγ , ulδ ) is the displacement of atom i (j, k, l) in the direction α (β , γ , δ ) , and ϕ , ψ and χ are the second-order harmonic, third-order and fourth-order anhamonic interatomic force constants, i.e. the second-, third- and fourth-order derivatives of the total energy with respect to the atom displacements from the equilibrium position. The unknown force constants can then be obtained by differentiating the displacement-force relation [53]:

Fiα = − ∑ ∑ ϕijαβ uiβ − i

β

1 2

∑ ∑ ψijkαβγ ujβ ukγ − jk

βγ

1 6

αβγδ β γ δ uj uk ul ∑ ∑ ψijkl jkl βγδ

(15) The interlayer coupling within the TiS2 is described by the LennardJones potential in Eq. (11). The empirical parameters for interlayer coupling are obtained by matching the potential surface as a function of interlayer distance. The optimized parameters are σTi − Ti = 2.8286 Å, ϵTi − Ti= 0.000737 eV, σS − S = 3.135 Å and ϵS − S= 0.0238 eV. However, MD simulations based on the full harmonic, third- and fourth-order force constant is too costly. To make the computation affordable, third

Fig. 11. Temperature-dependent thermal conductivity of MAPbI3 by MD simulations compared with the experimental results by A. Pisoni et al. [117] and direct AIMD simulation by Yue et al. [124].

Fig. 12. (a) phonon dispersion and partial density of states in pseudocubic MAPbI3. (b) frequency-dependent phonon group velocity. (c) Phonon lifetime spectrum as a function of frequency. (Reproduced with permission from Yue et al. [124], Copyright 2016 by American Physical Society).

403

Nano Energy 41 (2017) 394–407

X. Qian et al.

Fig. 13. (a) Thermal conductivity of pristine TiS2 and the organic intercalated hybrid TiS2 superlattice calculated by Green-Kubo integration (b–c). SED of the pristine and organic intercalated TiS2 superlattice. (Reproduced with permission from Wan et al. [13] Copyright 2015 by Nature Publishing Group).

interaction are also considered, respectively.

and fourth order constants for only two-body interactions are kept, and those of many body interactions are neglected [12]. For the organic part, CVFF potential [82] is adopted to describe the vibrations of the organic hexylammonium ions [CH3(CH2)5NH3]+. Due to the polar nature of the organic intercalant, the organic-inorganic couplings are modeled using both Lennard-Jones and Coulombic potentials which accounts both the van der Waals and the electrostatic interactions, respectively. When the organic ions are intercalated, the Ti atoms are reduced during the electrochemical process and one unit of negative charge is added to the Ti atom which is closest to the N atoms in the organic ions. The parameters for van der Waals coupling between the organic ions and the inorganic TiS2 can be obtained using the geometric mixing rules [83]. Fig. 13a shows the Green-Kubo integration of thermal conductivity for pristine TiS2 and the intercalated organicinorganic TiS2 superlattice. An order-of-magnitude reduction in the thermal conductivity is observed after the organic intercalation of TiS2. To explain the thermal conductivity reduction, SED analysis is conducted as shown in Fig. 13b–c. A broadening of the SED peaks is observed after TiS2 is intercalated with the organic ions, indicating stronger phonon scattering. The thermal conductivity reduction can be understood as a result of the inorganic TiS2 being ionically combined with the dangling hexylammonium ions. It is noted that the calculated in-plane thermal conductivity for both pristine TiS2 (13.2 W/m K) and the intercalated TiS2 (1.8 W/m K) is much higher than measurement (4.2 W/m K for TiS2 and 0.12 W/m K after organic intercalation) [12]. The reason causing the overestimated thermal conductivity could be attributed to two factors. One is that our calculations are for perfect single crystals while the experimental sample might inevitably contain some isotopes, vacancies, grain boundaries and other imperfectness, which might reduce the effective thermal conductivity considerably. In particular, for the intercalated TiS2, the wavy structure and local atomic disordering of the TiS2 layers were clearly identified from the STEM picture, which might explain the obtained ultralow thermal conductivity of intercalated TiS2 from the measurements. The other possibility is that the third- and fourth- order force constants corresponding to many body interactions are neglected in MD calculations to make computation affordable. Carefully including the components of the third- and fourth- force constants corresponding to many-body interaction into the force-field could substantially reduce the thermal conductivity. Phonon Boltzmann transport equation (BTE) calculation was conducted to predict the thermal conductivity of monolayer TiS2. The thermal conductivity of a 10 µm TiS2 monolayer is 24.7 W/m K when only two-body interaction, and it is reduced to 9.7 W/m K three-body

5. Conclusions and outlook In summary, based on the authors’ recent work, this article systematically presents an integrated first-principles-driven computational strategy for modeling the thermal conductivity of emerging hybrid organic-inorganic crystals and superlattices. This novel simulation strategy avoids both the computational challenges in direct first-principles simulation and the limits of classical molecular dynamics simulation. By either fitting the energy surface or interatomic forces from the ab-initio based first-principles simulations, empirical potential fields have been successfully developed for a range of emerging hybrid organic-inorganic crystals and superlattices, including organic-inorganic II–VI based hybrid crystals, the organometal halide perovskites, and the organic-intercalated TiS2 superlattice. With the potential fields developed, the classical molecular dynamics simulation is then used to predict the anisotropic thermal conductivity of these novel materials which are of great technical interests including photovoltaics, thermoelectrics and light emitting diodes. When implementing the DFT-MD method described in this Review, the readers should be aware that establishing the potential fields by matching ab-initio potential surface or inter-atomic forces are materialspecific and non-transferable. The potential field should be re-constructed when the atomic structure is changed, even with the same atomic species. For example, the potential fields need to be constructed separately for the tetragonal and pseudocubic phases of CH3NH3PbI3. Recently, the machine learning techniques based on Gaussian regression [125,126] and neural network [127] emerged as new methods for constructing potential fields with high accuracy and high transferability, which might be promising for modeling thermal transport in hybrid organic-inorganic materials in the near future. However, the machine learning-based potential is still limited to simple systems with no more than three elements, and there still lacks good method to describe long-range interactions [128], which are important in hybrid materials. This review provides a guidance for modeling thermal conductivity of materials with complex atomic structures from the first principles, which could be a good tutorial for both graduate students and experienced researchers who are interested in thermal energy transport in emerging hybrid materials.

404

Nano Energy 41 (2017) 394–407

X. Qian et al.

Acknowledgement This work was supported by NSF (Grant nos. 0846561 and 1512776). The DFT calculations are performed using the Summit Supercomputer, which is supported by NSF (awards ACI-1532235 and ACI-1532236), University of Colorado Boulder and Colorado State University. MD simulations are performed using the Oak Ridge Leadership Computing Facility operated by Oak Ridge National Laboratory.

[26]

[27]

[28]

[29]

Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.nanoen.2017.09.047.

[30]

[31]

References [32]

[1] L. Niu, et al., Controlled synthesis of organic/inorganic Van Der Waals solid for tunable light-matter interactions, Adv. Mater. 27 (2015) 7800–7808. [2] G. Kickelbick, Hybrid Materials: Synthesis, Characterization, and Applications, Wiley-VCH, Weinheim, Germany, 2007 (ISBN: 9783527610495). [3] X. Fang, M. Roushan, R. Zhang, J. Peng, H. Zeng, J. Li, Tuning and enhancing white light emission of Ii–Vi based inorganic–organic hybrid semiconductors as single-phased phosphors, Chem. Mater. 24 (2012) 1710–1717. [4] L. Dou, et al., Atomically thin two-dimensional organic-inorganic hybrid perovskites, Science 349 (2015) 1518–1521. [5] J.H. Im, C.R. Lee, J.W. Lee, S.W. Park, N.G. Park, 6.5% efficient perovskite quantum-dot-sensitized solar cell, Nanoscale 3 (2011) 4088–4093. [6] N.J. Jeon, J.H. Noh, Y.C. Kim, W.S. Yang, S. Ryu, S.I. Seok, Solvent engineering for high-performance inorganic-organic hybrid perovskite solar cells, Nat. Mater. 13 (2014) 897–903. [7] H.S. Kim, et al., Lead iodide perovskite sensitized all-solid-state submicron thin film mesoscopic solar cell with efficiency exceeding 9%, Sci. Rep. 2 (2012) 591. [8] J.H. Heo, et al., Efficient inorganic–organic hybrid heterojunction solar cells containing perovskite compound and polymeric hole conductors, Nat. Photonics 7 (2013) 486–491. [9] M. Liu, M.B. Johnston, H.J. Snaith, Efficient planar heterojunction perovskite solar cells by vapour deposition, Nature 501 (2013) 395–398. [10] S. Luo, W.A. Daoud, Recent progress in organic–inorganic halide perovskite solar cells: mechanisms and material design, J. Mater. Chem. A 3 (2015) 8992–9010. [11] J. Burschka, N. Pellet, S.J. Moon, R. Humphry-Baker, P. Gao, M.K. Nazeeruddin, M. Gratzel, Sequential deposition as a route to high-performance perovskite-sensitized solar cells, Nature 499 (2013) 316–319. [12] C. Wan, et al., Flexible N-type thermoelectric materials by organic intercalation of layered transition metal dichalcogenide Tis2, Nat. Mater. 14 (2015) 622–627. [13] C. Wan, Y. Kodama, M. Kondo, R. Sasai, X. Qian, X. Gu, K. Koga, K. Yabuki, R. Yang, K. Koumoto, Dielectric mismatch mediates carrier mobility in organicintercalated layered Tis2, Nano Lett. 15 (2015) 6302–6308. [14] Y. Zhou, F. Liu, H. Wang, Novel organic-inorganic composites with high thermal conductivity for electronic packaging applications: a key issue review, Polym. Compos. 38 (2017) 803–813. [15] G.M. Whitesides, J.P. Mathias, C.T. Seto, Molecular self-assembly and nanochemistry: a chemical strategy for the synthesis of nanostructures, Xcience 254 (1991) 1312–1319. [16] Y. Ni, S. Zheng, K. Nie, Morphology and thermal properties of inorganic–organic hybrids involving epoxy resin and polyhedral oligomeric silsesquioxanes, Polymer 45 (2004) 5557–5568. [17] Ordonez-Miranda, J., Ras, M.A., Wunderle, B., Volz, S., Modelling and Measurement of the Thermal conductivity of composites with silver particles. In: Proceeding of Therminic 2016 – 22nd International Workshop, Budapest/ Hungary, 2016, pp. 121–124. [18] J. Ordonez-Miranda, et al., Measurement and modeling of the effective thermal conductivity of sintered silver pastes, Int. J. Therm. Sci. 108 (2016) 185–194. [19] A.M. Marconnet, M.A. Panzer, K.E. Goodson, Thermal conduction phenomena in carbon nanotubes and related nanostructured materials, Rev. Mod. Phys. 85 (2013) 1295–1326. [20] A.M. Marconnet, N. Yamamoto, M.A. Panzer, B.L. Wardle, K.E. Goodson, Thermal conduction in aligned carbon nanotube–polymer nanocomposites with high packing density, ACS Nano 5 (2011) 4818–4825. [21] T. Luo, J.R. Lloyd, Non-equilibrium molecular dynamics study of thermal energy transport in Au–Sam–Au junctions, Int. J. Heat Mass Transf. 53 (2010) 1–11. [22] M.D. Losego, I.P. Blitz, R.A. Vaia, D.G. Cahill, P.V. Braun, Ultralow thermal conductivity in organoclay nanolaminates synthesized via simple self-assembly, Nano Lett. 13 (2013) 2215–2219. [23] M.D. Losego, M.E. Grady, N.R. Sottos, D.G. Cahill, P.V. Braun, Effects of chemical bonding on heat transport across interfaces, Nat. Mater. 11 (2012) 502–506. [24] L. Hu, L. Zhang, M. Hu, J.-S. Wang, B. Li, P. Keblinski, Phonon interference at selfassembled monolayer interfaces: molecular dynamics simulations, Phys. Rev. B (2010) 81. [25] S. Majumdar, J.A. Sierra-Suarez, S.N. Schiffres, W.L. Ong, C.F. Higgs III,

[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]

405

A.J. McGaughey, J.A. Malen, Vibrational mismatch of metal leads controls thermal conductance of self-assembled monolayer junctions, Nano Lett. 15 (2015) 2985–2991. W.L. Ong, S.M. Rupich, D.V. Talapin, A.J. McGaughey, J.A. Malen, Surface chemistry mediates thermal transport in three-dimensional nanocrystal arrays, Nat. Mater. 12 (2013) 410–415. W.-L. Ong, S. Majumdar, J.A. Malen, A.J.H. McGaughey, Coupling of organic and inorganic vibrational states and their thermal transport in nanocrystal arrays, J. Phys. Chem. C 118 (2014) 7288–7295. K.T. Regner, S. Majumdar, J.A. Malen, Instrumentation of broadband frequency domain thermoreflectance for measuring thermal conductivity accumulation functions, Rev. Sci. Instrum. 84 (2013) 064901. W.L. Ong, E.S. O'Brien, P.S. Dougherty, D.W. Paley, C. Fred Higgs III, A.J. McGaughey, J.A. Malen, X. Roy, Orientational order controls crystalline and amorphous thermal transport in superatomic crystals, Nat. Mater. 16 (2017) 83–88. X. Huang, M. Roushan, T.J. Emge, W. Bi, S. Thiagarajan, J.H. Cheng, R. Yang, J. Li, Flexible hybrid semiconductors with low thermal conductivity: the role of organic diamines, Angew. Chem. Int. Ed. Engl. 48 (2009) 7871–7874. X. Huang, J. Li, H. Fu, The first covalent organic-inorganic networks of hybrid chalcogenides: structures that may lead to a new type of quantum wells, J. Am. Chem. Soc. 122 (2000) 8789–8790. F. Brivio, A.B. Walker, A. Walsh, Structural and electronic properties of hybrid perovskites for high-efficiency thin-film photovoltaics from first-principles, arXiv:1309.4215v1 (2013). J. Feng, Mechanical properties of hybrid organic-inorganic Ch3nh3bx3 (B = Sn, Pb; X = Br, I) perovskites for solar cell absorbers, Apl. Mater. 2 (2014) 081801. J.M. Frost, K.T. Butler, F. Brivio, C.H. Hendon, M. van Schilfgaarde, A. Walsh, Atomistic origins of high-performance in hybrid halide perovskite solar cells, Nano Lett. 14 (2014) 2584–2590. Y. He, G. Galli, Perovskites for solar thermoelectric applications: a first principle study of Ch3nh3ai3 (a = Pb and Sn), Chem. Mater. 26 (2014) 5394–5400. J. Liu, B. Yoon, E. Kuhlmann, M. Tian, J. Zhu, S.M. George, Y.C. Lee, R. Yang, Ultralow thermal conductivity of atomic/molecular layer-deposited hybrid organic-inorganic zincone thin films, Nano Lett. 13 (2013) 5594–5599. J. Liu, J. Zhu, M. Tian, X. Gu, A. Schmidt, R. Yang, Simultaneous measurement of thermal conductivity and heat capacity of bulk and thin film materials using frequency-dependent transient thermoreflectance method, Rev. Sci. Instrum. 84 (2013) 034902. X. Gu, R. Yang, First-principles prediction of phononic thermal conductivity of silicene: a comparison with graphene, J. Appl. Phys. 117 (2015) 025102. X. Gu, R. Yang, Phonon transport in single-layer transition metal dichalcogenides: a first-principles study, Appl. Phys. Lett. 105 (2014) 131903. X. Li, K. Maute, M.L. Dunn, R. Yang, Strain EFfects on the Thermal Conductivity of Nanostructures, Phys. Rev. B (2010) 81. X. Li, R. Yang, Equilibrium molecular dynamics simulations for the thermal conductivity of Si/Ge nanocomposites, J. Appl. Phys. 113 (2013) 104306. J. Liu, R. Yang, Tuning the thermal conductivity of polymers with mechanical strains, Phys. Rev. B (2010) 81. X. Qian, X. Gu, M.S. Dresselhaus, R. Yang, Anisotropic tuning of graphite thermal conductivity by lithium intercalation, J. Phys. Chem. Lett. 7 (2016) 4744–4750. D.G. Cahill, et al., Nanoscale thermal transport. II. 2003–2012, Appl. Phys. Rev. 1 (2014) 011305. D.G. Cahill, W.K. Ford, K.E. Goodson, G.D. Mahan, A. Majumdar, H.J. Maris, R. Merlin, S.R. Phillpot, Nanoscale thermal transport, J. Appl. Phys. 93 (2003) 793. T. Luo, G. Chen, Nanoscale heat transfer–from computation to experiment, Phys. Chem. Chem. Phys. 15 (2013) 3389–3412. G. Chen, Nanoscale Energy Transport and Conversion, Oxford University Press, 2005. X. Gu, R. Yang, Phonon transport and thermal conductivity in two-dimensional materials, Annu. Rev. Heat Transf. 19 (2016) 1–65. A.J.H. McGaughey, J.M. Larkin, Predicting phonon properties from equilibrium molecular dynamics simulations, Annu. Rev. Heat Transf. 17 (2014) 49–87. D.A. Broido, M. Malorny, G. Birner, N. Mingo, D.A. Stewart, Intrinsic Lattice thermal conductivity of semiconductors from first principles, Appl. Phys. Lett. 91 (2007) 231922. L. Lindsay, D.A. Broido, T.L. Reinecke, Thermal conductivity and large isotope effect in gan from first principles, Phys. Rev. Lett. 109 (2012) 095901. A. Ward, D.A. Broido, Intrinsic phonon relaxation times from first-principles studies of the thermal conductivities of Si and Ge, Phys. Rev. B (2010) 81. K. Esfarjani, H.T. Stokes, Method to extract anharmonic force constants from first principles calculations, Phys. Rev. B (2008) 77. W. Li, N. Mingo, L. Lindsay, D.A. Broido, D.A. Stewart, N.A. Katcho, Thermal conductivity of diamond nanowires from first principles, Phys. Rev. B (2012) 85. A. Ward, D.A. Broido, D.A. Stewart, G. Deinzer, Ab initio theory of the lattice thermal conductivity in diamond, Phys. Rev. B (2009) 80. K. Esfarjani, G. Chen, H.T. Stokes, Heat transport in silicon from first-principles calculations, Phys. Rev. B (2011) 84. G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, F. Mauri, Thermal conductivity of graphene and graphite: collective excitations and mean free paths, Nano Lett. 14 (2014) 6109–6114. L. Lindsay, D.A. Broido, N. Mingo, Flexural phonons and thermal transport in graphene, Phys. Rev. B (2010) 82. L. Lindsay, D.A. Broido, N. Mingo, Flexural phonons and thermal transport in multilayer graphene and graphite, Phys. Rev. B (2011) 83.

Nano Energy 41 (2017) 394–407

X. Qian et al. [60] H. Xie, M. Hu, H. Bao, Thermal conductivity of silicene from first-principles, Appl. Phys. Lett. 104 (2014) 131906. [61] A. Jain, A.J. McGaughey, Strongly anisotropic in-plane thermal transport in singlelayer black phosphorene, Sci. Rep. 5 (2015) 8501. [62] G. Qin, Q.B. Yan, Z. Qin, S.Y. Yue, M. Hu, G. Su, Anisotropic intrinsic lattice thermal conductivity of phosphorene from first principles, Phys. Chem. Chem. Phys. 17 (2015) 4854–4858. [63] B. Sun, X. Gu, Q. Zeng, X. Huang, Y. Yan, Z. Liu, R. Yang, Y.K. Koh, Temperature dependence of anisotropic thermal-conductivity tensor of bulk black phosphorus, Adv. Mater. (2017) 29. [64] J. Zhu, et al., Revealing the origins of 3d anisotropic thermal conductivities of black phosphorus, Adv. Electron. Mater. 2 (2016) 1600040. [65] X. Gu, B. Li, R. Yang, Layer thickness-dependent phonon properties and thermal conductivity of Mos2, J. Appl. Phys. 119 (2016) 085106. [66] X. Gu, R. Yang, Phonon transport in single-layer Mo1-xWxS2 alloy embedded with WS2 nanodomains, Phys. Rev. B (2016) 94. [67] L. Lindsay, D.A. Broido, T.L. Reinecke, First-principles determination of ultrahigh thermal conductivity of boron arsenide: a competitor for diamond? Phys. Rev. Lett. 111 (2013) 025901. [68] A. Mattoni, A. Filippetti, M.I. Saba, P. Delugas, Methylammonium rotational dynamics in lead halide perovskite by classical molecular dynamics: the role of temperature, J. Phys. Chem. C 119 (2015) 17421–17428. [69] T. Chen, B.J. Foley, B. Ipek, M. Tyagi, J.R. Copley, C.M. Brown, J.J. Choi, S.H. Lee, Rotational dynamics of organic cations in the Ch3nh3pbi3 perovskite, Phys. Chem. Chem. Phys. 17 (2015) 31278–31286. [70] A.M. Leguy, et al., The dynamics of methylammonium ions in hybrid organicinorganic perovskite solar cells, Nat. Commun. 6 (2015) 7124. [71] S. Stackhouse, L. Stixrude, B.B. Karki, Thermal conductivity of periclase (Mgo) from first principles, Phys. Rev. Lett. 104 (2010) 208501. [72] J. Haile, Molecular Dynamics Simulation 18 Wiley, New York, 1992. [73] P.C. Howell, Comparison of molecular dynamics methods and interatomic potentials for calculating the thermal conductivity of silicon, J. Chem. Phys. 137 (2012) 224111. [74] M.B. Kanoun, A.E. Merad, H. Aourag, J. Cibert, G. Merad, Molecular-dynamics simulations of structural and thermodynamic properties of znte using a three-body potential, Solid State Sci. 5 (2003) 1211. [75] L. Lindsay, D.A. Broido, Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene, Phys. Rev. B (2010) 81. [76] X.W. Zhou, D.K. Ward, J.E. Martin, F.B. van Swol, J.L. Cruz-Campa, D. Zubia, Stillinger-Weber potential for the II–VI elements Zn-Cd-Hg-S-Se-Te, Phys. Rev. B (2013) 88. [77] D.K. Ward, X.W. Zhou, B.M. Wong, F.P. Doty, J.A. Zimmerman, Analytical bondorder potential for the cadmium telluride binary system, Phys. Rev. B (2012) 85. [78] D.K. Ward, X.W. Zhou, B.M. Wong, F.P. Doty, J.A. Zimmerman, Analytical bondorder potential for the Cd-Zn-Te ternary system, Phys. Rev. B (2012) 86. [79] Y.T. Cheng, T. Liang, J.A. Martinez, S.R. Phillpot, S.B. Sinnott, A. Charge, Optimized many-body potential for titanium nitride (Tin), J. Phys. Condens. Matter.: Inst. Phys. J. 26 (2014) 265004. [80] Y.T. Cheng, T.R. Shan, T. Liang, R.K. Behera, S.R. Phillpot, S.B. Sinnott, A. Charge, Optimized many-body (comb) potential for titanium and titania, J. Phys. Condens. Matter.: Inst. Phys. J. 26 (2014) 315007. [81] X. Qian, X. Gu, R. Yang, Anisotropic thermal transport in organic–inorganic hybrid crystal Β-Znte(En)0.5, J. Phys. Chem. C (2015). [82] J.-R. Hill, J. Sauer, Molecular mechanics potential for silica and zeolite catalysits based on ab initio calculations. 2. Aluminosilicates, J. Phys. Chem. 99 (1995) 9536–9550. [83] H. Sun, Compass: an ab initio force-field optimized for condensed-phase applications-overview with details on alkane and benzene compounds, J. Phys. Chem. B 102 (1998) 7338–7364. [84] P. Dauber-Osguthorpe, V.A. Roberts, D.J. Osguthorpe, J. Wolff, M. Genest, A.T. Hagler, Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system, Proteins 4 (1988) 31–47. [85] J.E.D. Jr, J.J. More, Quasi-newton methods, motivation and theory, SIAM Rev. (1977) 19. [86] Rohskopf, A., Seyf, H.R., Gordiz, K., AsegunHenry, Phonon_Optimized_Potential. Pdf. arXiv:1610.02353, 2017. [87] F. Müller-Plathe, A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity, J. Chem. Phys. 106 (1997) 6082–6085. [88] D.P. Sellan, E.S. Landry, J.E. Turney, A.J.H. McGaughey, C.H. Amon, Size effects in molecular dynamics thermal conductivity predictions, Phys. Rev. B (2010) 81. [89] E. Lampin, P.L. Palla, P.A. Francioso, F. Cleri, Thermal conductivity from approach-to-equilibrium molecular dynamics, J. Appl. Phys. 114 (2013) 033525. [90] P.K. Schelling, S.R. Phillpot, P. Keblinski, Comparison of atomic-level simulation methods for computing thermal conductivity, Phys. Rev. B (2002) 65. [91] A.J.H. McGaughey, Phonon Transport in Molecular Dynamics Simulations: Formulation and Thermal Conductivity Prediction, University of Michigan, 2004. [92] W.J. Parker, R.J. Jenkins, C.P. Butler, G.L. Abbott, Flash method of determining thermal diffusivity, heat capacity, and thermal conductivity, J. Appl. Phys. 32 (1961) 1679–1684. [93] M. Kaviany, Heat Transfer Physics, Cambridge University Press, 2008. [94] A.J.H. McGaughey, M. Kaviany, Phonon transport in molecular dynamics simulations: formulation and thermal conductivity prediction, Adv. Heat Transf. 39 (2006) 169–255. [95] Z. Wang, S. Safarkhani, G. Lin, X. Ruan, Uncertainty quantification of thermal

[96]

[97] [98]

[99] [100] [101] [102]

[103] [104]

[105] [106] [107]

[108] [109] [110]

[111]

[112]

[113] [114]

[115]

[116]

[117]

[118] [119] [120]

[121]

[122]

[123]

[124]

[125] [126]

406

conductivities from equilibrium molecular dynamics simulations, Int. J. Heat Mass Transf. 112 (2017) 267–278. J. Turney, E. Landry, A. McGaughey, C. Amon, Predicting phonon properties and thermal conductivity from anharmonic lattice dynamics calculations and molecular dynamics simulations, Phys. Rev. B (2009) 79. T. Feng, B. Qiu, X. Ruan, Anharmonicity and necessity of phonon eigenvectors in the phonon normal mode analysis, J. Appl. Phys. 117 (2015) 195102. J.A. Thomas, J.E. Turney, R.M. Iutzi, C.H. Amon, A.J.H. McGaughey, Predicting phonon dispersion relations and lifetimes from the spectral energy density, Phys. Rev. B (2010) 81. X. Qian, X. Gu, R. Yang, Lattice thermal conductivity of organic-inorganic hybrid perovskite Ch3nh3pbi3, Appl. Phys. Lett. 108 (2016) 063902. W. Lv, A. Henry, Direct calculation of modal contributions to thermal conductivity via Green–Kubo modal analysis, New J. Phys. 18 (2016) 013028. M.T. Dove, Introduction to Lattice Dynamics 4 Cambridge University Press, 1993. B.L. Huang, A.J.H. McGaughey, M. Kaviany, Thermal conductivity of metal-organic framework 5 (Mof-5): Part I. Molecular dynamics simulations, Int. J. Heat Mass Transf. 50 (2007) 393–404. B. Huang, M. Kaviany, Ab initio and molecular dynamics predictions for electron and phonon transport in bismuth telluride, Phys. Rev. B 77 (2008) 125209. B. Huang, M. Kaviany, Filler-reduced phonon conductivity of thermoelectric skutterudites: ab initio calculations and molecular dynamics simulations, Acta Mater. 58 (2010) 4516–4526. J.E. Jones, On the determination of molecular fields. II. From the equation of state of a gas, Proc. R. Soc. A 106 (1924) 463. D. Connétable, O. Thomas, First-principles study of the structural, electronic, vibrational, and elastic properties of orthorhombic NiSi, Phys. Rev. B (2009) 79. C.-Y. Moon, G.M. Dalpian, Y. Zhang, S.-H. Wei, X.-Y. Huang, J. Li, Study of phase selectivity of organic−inorganic hybrid semiconductors, Chem. Mater. 18 (2006) 2805–2809. Y. Zhang, et al., Zero thermal expansion in a nanostructured inorganic-organic hybrid crystal, Phys. Rev. Lett. (2007) 99. J.W. Jaeken, S. Cottenier, Solving the christoffel equation: phase and group velocities, Comput. Phys. Commun. 207 (2016) 445–451. C.C. Stoumpos, C.D. Malliakas, M.G. Kanatzidis, Semiconducting tin and lead iodide perovskites with organic cations: phase transitions, high mobilities, and nearinfrared photoluminescent properties, Inorg. Chem. 52 (2013) 9019–9038. A. Kojima, K. Teshima, Y. Shirai, T. Miyasaka, Organnometal halide perovskites as visibile-light sensitizers for photovoltaic cells, J. Am. Chem. Soc. 131 (2009) 6050–6051. T. Katayama, K. Setoura, D. Werner, H. Miyasaka, S. Hashimoto, Picosecond-tonanosecond dynamics of plasmonic nanobubbles from pump-probe spectral measurements of aqueous colloidal gold nanoparticles, Langmuir: ACS J. Surf. Colloids 30 (2014) 9504–9513. D. Shi, et al., Low trap-state density and long carrier diffusion in organolead trihalide perovskite single crystals, Science (2015) 347. S.D. Stranks, G.E. Eperon, G. Grancini, C. Menelaou, M.J.P. Alcocer, Tomas Leijtens, L.M. Herz, A. Petrozza, H.J. Snaith, Electron-hole diffusion lengths exceeding 1 μm in an organometal trihalide perovskite absorber, Science (2013) 342. C. Lee, J. Hong, A. Stroppa, M.-H. Whangbo, J.H. Shim, Organic–inorganic hybrid perovskites Abi3(a = Ch3nh3, Nh2chnh2; B = Sn, Pb) as potential thermoelectric materials: a density functional evaluation, RSC Adv. 5 (2015) 78701–78707. X. Mettan, R. Pisoni, P. Matus, A. Pisoni, J. Jaćimović, B. Náfrádi, M. Spina, D. Pavuna, L. Forró, E. Horváth, Tuning of the thermoelectric figure of merit of Ch3nh3mi3(M = Pb,Sn) photovoltaic perovskites, J. Phys. Chem. C 119 (2015) 11506–11510. A. Pisoni, J. Jaćimović, O.S. Barišić, M. Spina, R. Gaál, L. Forró, E. Horváth, Ultralow thermal conductivity in organic–inorganic hybrid perovskite Ch3nh3pbi3, J. Phys. Chem. Lett. 5 (2014) 2488–2492. R.A. Buckingham, The classical equation of state of gaseous helium, neon and argon, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 168 (1938) 264–283. W. Tang, E. Sanville, G. Henkelman, A. Grid-Based Bader, Analysis algorithm without lattice bias, J. Phys. Condens. Matter.: Inst. Phys. J. 21 (2009) 084204. F. Brivio, J.M. Frost, J.M. Skelton, A.J. Jackson, O.J. Weber, M.T. Weller, A.R. Goñi, A.M.A. Leguy, P.R.F. Barnes, A. Walsh, Lattice dynamics and vibrational spectra of the orthorhombic, tetragonal, and cubic phases of methylammonium lead iodide, Phys. Rev. B (2015) 92. T. Hata, G. Giorgi, K. Yamashita, The effects of the organic-inorganic interactions on the thermal transport properties of Ch3nh3pbi3, Nano Lett. 16 (2016) 2749–2753. M. Wang, S. Lin, Anisotropic and ultralow phonon thermal transport in organicinorganic hybrid perovskites: atomistic insights into solar cell thermal management and thermoelectric energy conversion efficiency, Adv. Funct. Mater. 26 (2016) 5297–5306. C. Caddeo, C. Melis, M.I. Saba, A. Filippetti, L. Colombo, A. Mattoni, Tuning the thermal conductivity of methylammonium lead halide by the molecular substructure, Phys. Chem. Chem. Phys. 18 (2016) 24318–24324. S.-Y. Yue, X. Zhang, G. Qin, J. Yang, M. Hu, Insight into the collective vibrational modes driving ultralow thermal conductivity of perovskite solar cells, Phys. Rev. B (2016) 94. A.P. Bartók, G. Csányi, Gaussian approximation potentials: a brief tutorial introduction, Int. J. Quantum Chem. 115 (2015) 1051–1057. A.P. Bartok, M.C. Payne, R. Kondor, G. Csanyi, Gaussian approximation potentials: the accuracy of quantum mechanics, without the electrons, Phys. Rev. Lett. 104 (2010) 136403.

Nano Energy 41 (2017) 394–407

X. Qian et al.

Dr. Ronggui Yang (Ph.D. MIT 2006, M.S. UCLA 2001, M.S. Tsinghua University 1999, B.S. Xi’an Jiaotong University 1996) is a Professor of Mechanical Engineering at the University of Colorado Boulder. His research interests are on the fundamentals of micro/ nanoscale transport phenomena and the applications of micro/nano technologies for energy conversion, storage and thermal management. He is an ASME Fellow. He currently chairs the ASME Heat Transfer Division's K-9 Technical Committee on Nanoscale Transport Phenomena and serves as an Associate Editor of ASME Journal of Heat Transfer, among other service activities.

[127] N. Artrth, J. Behler, High-dimensional network potentials for metal surfaces: a prototype study for copper, Phys. Rev. B 85 (2012) 045439. [128] A.P. Bartók, R. Kondor, G. Csányi, On representing chemical environments, Phys. Rev. B (2013) 87.

Xin Qian is currently a Ph.D. student in Mechanical Engineering at University of Colorado Boulder with Professor Ronggui Yang. He received his B.S. at Huazhong University of Science and Technology in 2014. His research focused on modeling and experimental characterization of thermal transport in complex materials including hybrid organic-inorganic materials, high temperature thermal barrier materials and wide-bandgap semiconductors.

Xiaokun Gu is currently an assistant professor of Institute of Engineering Thermophysics, Shanghai Jiao Tong University. He received his Ph.D. degree in 2015 at University of Colorado Boulder with Professor Ronggui Yang. His research focused on theoretical modeling of phonon transport in hybrid organic-inorganic materials and low dimensional materials.

407