PRL 112, 065901 (2014)
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
Lattice Thermal Conductivity of Si1−x Gex Nanocomposites Claudio Melis and Luciano Colombo* Dipartimento di Fisica, Università di Cagliari, Cittadella Universitaria, I-09042 Monserrato (Ca), Italy (Received 24 July 2013; published 11 February 2014) We calculate the lattice thermal conductivity in model Si1−x Gex nanocomposites by molecular dynamics in a transient thermal conduction regime. Our simulations provide evidence that thermal transport depends only marginally on stoichiometry in the range 0.2 ≤ x ≤ 0.8, while it is deeply affected by the granulometry. In particular, we show that Si1−x Gex nanocomposites have lattice thermal conductivity below the corresponding bulk alloy with the same stoichiometry. The main role in affecting thermal conduction is provided by grain boundaries, which largely affect vibrational modes with a long mean-free path. DOI: 10.1103/PhysRevLett.112.065901
PACS numbers: 66.70.-f, 81.07.Bc, 44.10.+i
The physics of thermoelectric conversion is summarized by the figure of merit ZT ¼ S2 Tσ=κ characterizing the efficiency of heat conversion to electricity at temperature T (here, S is the Seebeck coefficient, while σ and κ are the electrical and thermal conductivity, respectively). Energy applications call for high ZT values or, equivalently, for minimal κ (still preserving good charge transport characteristics) [1]. In most systems of current interest, the phonon (lattice) contribution to κ is largely dominant to the electron one. This motivates the search for new materials with very low lattice thermal conductivity, a situation obtained either by tailoring new chemically complex materials with suitable phonon properties [2] or by nanostructuring more conventional materials so that internal boundaries (phase or grain boundaries, as well as interfaces) maximize phonon scattering [3,4]. In this Letter we offer a contribution to this ongoing effort by discussing, in particular, the lattice thermal conduction in nanocrystalline Si1−x Gex nanocomposites (nc-Si1−x Gex ). To this aim, we exploit a procedure to calculate κ through a MD simulation in a transient conduction (nonequilibrium) condition, which has appealing features with respect to standard techniques. By establishing a constant thermal gradient ∂T=∂z within an homogeneous and incompressible material, the corresponding steady-state thermal conduction regime is well described by the Fourier equation jz ¼ −κ∂T=∂z, where jz is the heat current [5]. By MD simulations κ is usually computed either by the Green-Kubo formalism [6] or by so-called direct methods [7,8]. In the first case, the statistical theory of transport coefficients is used which provides κ through the R current-current autocorrelation function κ ¼ ð1=kB T 2 VÞ 0þ∞ hjz ðtÞjz ð0Þidt, where kB is the Boltzmann constant, V is the system volume, and h⋅i indicates the ensemble average. Direct methods, instead, compute separately the external perturbation ∂T=∂z and the system response jz and straightforwardly provide the thermal conductivity as the response-to-perturbation ratio κ ¼ −hjz ðtÞi=∂T=∂z. Practical applications of both 0031-9007=14=112(6)=065901(5)
approaches are limited by numerical as well as conceptual restrictions. First of all, we remark that the heat current for a system of N interacting atoms is typically written as P ⃗ jðtÞ ¼ dtd Ni¼1 ϵi ðtÞ⃗ri ðtÞ, where r⃗ i and ϵi are the position and energy density of the ith particle [8,9]. Such an expression is very cumbersome to calculate except for the simplest interatomic potentials and is not directly applicable when using total-energy functionals. In addition, both hjz ðtÞjz ð0Þi and hjz ðtÞi are very much slowly convergent in time and/or size [10]. Direct methods are additionally affected by the very long simulation times needed to establish a steady-state condition, once a given temperature gradient or heat current is imposed. This often makes such methods impractical, especially for systems (as those here investigated) that, by their very structure, require large simulation cells and/or are characterized by low thermal conduction. In this work we follow a different MD approach, which simulates a transient thermal conduction regime [11,12] providing κ by directly solving the heat equation ∂T ∂ 2T ¼ κ¯ 2 ; ∂t ∂z
(1)
where κ¯ ¼ κ=ρCv is the thermal diffusivity of the system with density ρ and specific heat Cv . This method is referred to as approach-to-equilibrium MD (AEMD) [12] and it is adopted here to investigate thermal conduction in ncSi1−x Gex systems differing by stoichiometry and granulometry. We provide evidence that nc-Si1−x Gex systems have room-temperature thermal conductivity below the bulk alloy limit [13]. Furthermore, we show that κ is largely affected by granulometry, while it only marginally depends on stoichiometry over an extended range of x values. Border conditions must be imposed on Eq. (1), specifying the initial temperature profile Tðz; t ¼ 0Þ along the z direction of a simulation cell as large as Lz . Within a MD simulation the natural choice is selecting periodic
065901-1
© 2014 American Physical Society
conditions and, therefore, we assume a steplike periodic profile with Tðz; t ¼ 0Þ ¼ T 1 in the semicell 0 ≤ z < Lz =2 and Tðz; t ¼ 0Þ ¼ T 2 in the semicell Lz =2 ≤ z < Lz (with T 1 > T 2 ). It easy to prove that the general solution of Eq. (1) can be cast in the form Tðz; tÞ ¼
∞ T1 þ T2 X 2 Bn sinðαn zÞe−αn κ¯ t ; þ 2 n¼1
(2)
where αn ¼ 2πn=Lz and Bn are the coefficients of the Fourier expansion of the periodic Tðz; tÞ function. By aging the system in a microcanonical MD simulation, the initial steplike temperature profile is progressively smoothed by thermal conduction; the average temperatures hT 1 i and hT 2 i in the two semicells vary in time and the systems eventually approach a uniform temperature (i.e., equilibrium) configuration. During such a transient regime we can define the time-dependent difference in average temperatures ΔTðtÞ ¼ hT 1 ðtÞi − hT 2 ðtÞi which, through Eq. (2), is easily calculated analytically as ΔTðtÞ ¼
∞ X
2
Cn e−αn κ¯ t ;
(3)
n¼1
where Cn ¼ 8ðT 1 − T 2 Þ½cosðαn Lz =2Þ − 12 =α2n L2z . This is the key equation of AEMD: simply, the ΔTðtÞ profile must be computed on the fly during the microcanonical run and then fitted by Eq. (3) to get κ¯ . Quantum corrections to the specific heat in the corresponding thermal conductivity κ ¼ ρCv κ¯ must be duly considered below the Debye temperature of the investigated system [14]. While more details of the method are discussed elsewhere [12], we note that AEMD provides evidence that Eq. (1) is still valid at the length scales typically explored by MD simulations. Figure 1 shows the time evolution of a steplike temperature profile with ΔTð0Þ ¼ 200 K generated across a Si0.5 Ge0.5 bulk alloy sample with Lz ¼ 166.35 nm. While the smooth lines represent the 500
formal solutions of Eq. (1), the noisy lines represent the direct MD calculation of local temperature during the microcanonical run. We point out other important AEMD features, namely, as anticipated, there is no need to calculate or impose the heat current or to run long simulations to establish a steady state. The method is also very robust: no result discussed in this Letter is significantly affected by varying the initial temperature step ΔTð0Þ in the range 50–200 K or by varying the number of exponential terms in Eq. (2) in the range 10–100. In short, AEMD is characterized by a clean theoretical foundation, ease of implementation, numerical robustness, and comparatively light computational effort. By using AEMD with the Tersoff force field [15,16] (which provides good predictions on the structure of Si1−x Gex alloys [17,18], as well as reasonable results for their thermal properties [10,12]), we first calculated the room-temperature thermal conductivity of Si1−x Gex bulk alloys. In order to predict the actual values for bulk samples, for any chemical composition of the alloy, we simulated several periodically repeated simulation cells with the same lateral dimensions, 3.8 nm ≤ Lx ¼ Ly ≤ 3.9 nm (varying with the stoichiometry), and z thickness varying in the range 273 ≤ Lz ≤ 563 nm. By the standard procedure [8] of plotting 1=κ versus 1=Lz , we obtained the bulk κ value by linearly extrapolating the calculated data in the limit 1=Lz → 0. The initial steplike temperature profile was set by Nosé-Hoover thermostatting with ΔTð0Þ ¼ 200 K, while in the following microcanonical MD run the equations of motion were integrated by the velocity Verlet algorithm with a 1 fs time step. The LAMMPS code was used [19]. Our results are compared in Fig. 2 with state-of-the-art ab initio data [13] and available experiments [20–22]. The resulting picture represents the reliability of the present method: κ is basically independent of the stoichiometry in the range 0.2 ≤ x ≤ 0.8, but its value sharply increases toward the c-Si and c-Ge values as x → 0 or 1, respectively. In particular, we obtained 216 18 and 65 6 W K−1 m−1 250
t=0 ps t=20 ps t=50 ps t=100 ps
450
This work Ref. [13] Ref. [20] Ref. [21] Ref. [22]
200
400
-1 -1 κ [W K m ]
Temperature [K]
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
PRL 112, 065901 (2014)
350 300 250
150 100 50
200 150
0 0
20
40
60
80 100 120 140 160
0
z [nm]
0.2
0.4
0.6
0.8
1
% Germanium Content
FIG. 1 (color online). Evolution of the temperature profile across a Si0.5 Ge0.5 bulk alloy sample with Lz ¼ 166.35 nm as provided by solving analytically Eq. (1) (smooth lines) and by direct MD calculation of local temperature (noisy color lines).
FIG. 2 (color online). Room-temperature thermal conductivity of bulk Si1−x Gex alloys predicted by AEMD simulations (red squares). Experimental data (color triangles and diamonds) and ab initio calculations (black dots) are shown for comparison.
065901-2
PRL 112, 065901 (2014)
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
for c-Si and c-Ge to compare with the corresponding range of experimental values 130–240 and 50–60 W K−1 m−1 , respectively (see Ref. [10] and references therein). While the agreement for Ge is very good, we remark that our result for Si is better than most of the previously reported values obtained with other methods, by using the same Tersoff force field. Only the calculation by He et al. [10] obtained a slightly better κ for c-Si, using, however, a huge simulation cell and extremely long simulation times. Generating atomistic models of nc-Si1−x Gex is a challenging and very computer-intensive problem, since we need to correctly reproduce both local-scale and global structural features (respectively, grain boundary atomic architecture and grain size or shape distribution) for many different stoichiometries. With this in mind, we followed a multistep procedure repeatedly applied to simulation cells that, by construction, initially contained just Si atoms and had lateral dimensions Lx ¼ 2.715 nm and Ly ¼ 27.15 nm. The z thickness was instead varied in the range 54.3 ≤ Lz ≤ 135.7 nm for the reason discussed above. Each system was at first fully amorphized by quenching from the melt. Then, N g sites were selected at random in the yz plane and around each of them a cylindrical void (passing across the full Lx dimension) was created by removing atoms. Voids were distributed avoiding overlap, and their initial radius was randomly set above the capillarity threshold [23] (∼0.5 nm). Next, each void was filled by a crystalline cylinder randomly rotated in the yz plane. The resulting mixed structures (an amorphous matrix spotted by crystalline seeds) were annealed at constant T ¼ 1200 K for very long time, varying from 1.5 to 4.0 ns for increasing granularity. During this step seed grains underwent growth until nanocrystalline silicon systems (hereafter referred to as templates) were eventually obtained. Three different templates are shown in Fig. 3: the resulting granularity hdg i reported in the figure for each sample
was effectively controlled by the initial choice of N g . At this stage nc-Si1−x Gex samples were obtained by randomly replacing Si atoms with Ge atoms in the templates, followed by further careful relaxation. This alchemy was flanked by a self-affine rescaling of the atomic positions, so that the actual lattice spacing of any nc-Si1−x Gex sample was set as aðSi1−x Gex Þ ¼ ð1 − xÞaðSiÞ þ xaðGeÞ, where aðSiÞ and aðGeÞ are the Si and Ge lattice constant provided by the Tersoff potential, respectively. This somewhat complex procedure is not only effective in producing trustworthy atomistic models (compare Fig. 3 with TEM pictures reported in Refs. [22,24]), but it is also very efficient: as a matter of fact, the most computerintensive step (i.e., the growth of crystalline grains against amorphous regions) must be followed just once, while Si → Ge replacements and following relaxation are comparatively light and, therefore, can be effortlessly repeated for any stoichiometry we are interested in. In total we generated a set of five different templates, with 3 ≤ hdg i ≤ 25 nm, corresponding to a number of atoms varying in the range ð2–5Þ × 105 (according to the actual Lz value of the simulation cell), indeed a large enough number to warrant that the random chemical arrangement of Si and Ge atoms is representative of real samples at any investigated stoichiometry. More importantly, however, the explored range of average grain size perfectly matches real nc-Si1−x Gex samples [22,24]. The values of thermal conductivity (for any given stoichiometry and granularity) discussed below were obtained as the 1=Lz → 0 extrapolation of a set of κ values corresponding to the five different systems previously generated. Furthermore, the thermal conductivity of each single system was calculated as the configurational average over two different samples with the same Lz , x, and hdg i. We first address the key question of whether thermal conductivity is larger in bulk alloys or in nanocomposites. In Fig. 4 the remarkable result of the present investigation is reported for a nc-Si1−x Gex sample of granularity hdg i ¼ 10.5 nm. Here we provide evidence that κ is smaller by a ∼50% in nanocomposites, consistent with experimental data which report a very similar reduction in thermal 6
Bulk alloy Nanocomposite
κ [W K -1 m-1]
5 4 3 2 1
FIG. 3 (color online). Left: Snapshots of the final configuration obtained after the growth process described in the text for three template samples differing for average grain size hdg i. Right: Enlargement of the marked region, showing the atomistic structure of grain boundaries.
0.2
0.3
0.4
0.5
0.6
0.7
0.8
% Germanium Content
FIG. 4 (color online). Predicted room-temperature thermal conductivity of nc-Si1−x Gex samples with hdg i ¼ 10.5 nm and different stoichiometry.
065901-3
90
4 Bulk alloy Nanocomposite
80
3
70 60 50 40 30
2.5 2 1.5 1
20
0.5
10 0 0.0001 0.001
Nanocomposite Bulk alloy
3.5
-1 κ [W K m-1]
Thermal conductivity accumulation [%]
100
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
PRL 112, 065901 (2014)
0.01
0.1
1
0
10
0
5
10
15
20
25
30
Average grain size [nm]
Lz [µm]
FIG. 5 (color online). Thermal conductivity accumulation function for a nc-Si0.5 Ge0.5 sample (black dots) and a Si0.5 Ge0.5 bulk alloy (red squares). Solid lines show the corresponding accumulation extracted from the linear 1=κ versus 1=Lz extrapolation.
conductivity of both p- and n-type samples [22,24]. While experiments were collected for the sole x ¼ 0.2 stoichiometry, Fig. 4 shows that this feature is indeed valid for any chemical composition. This is an important evidence supporting the use of nc-Si1−x Gex systems as thermoelectric materials of possibly large technological impact. The reduced value of κ is clearly due to the enhancement of scattering phenomena provided by the network of grain boundaries (see Fig. 3). This statement is consistent with recent results on short-period Si=Ge superlattices [25,26], where the major role played by interface scattering has been clearly defined. We further substantiate this concept by the calculation of the thermal conductivity accumulation function, here defined as the ratio κðLz Þ=κ ∞ between its value computed for a simulation cell with length Lz and its extrapolated value for a bulklike (infinite) sample. κðLz Þ=κ ∞ gives the contribution to the thermal conductivity provided by vibrational modes with mean-free path λ up to Lz . In Fig. 5 it is shown that for a nc-Si0.5 Ge0.5 sample, as much as ∼80% of the thermal conductivity is provided by vibrational modes with λ ≤ 100 nm. On the other hand, in the bulk alloy the same contribution is provided by vibrations with much longer λ, up to 1 μm, consistent with the estimation by Garg et al. [13]. Another feature emerging from Fig. 4 is that, similarly to bulk alloys, the actual chemical composition seems to play a minor role: as a matter of fact, κ is basically independent of stoichiometry, although there is an indication that a small Ge content should be preferred in order to fully minimize κ. We can then draw the conclusion that alloying is not the ultimate solution for minimizing thermal conduction in Si1−x Gex systems; rather, this role is played by grain boundaries. In order to better characterize the role of the nanostructure, we have extended the above calculations to the full set of nc-Si1−x Gex samples with x ¼ 0.5. Figure 6 shows that κ increases monotonically with granulometry, and, for a large enough average grain size, the thermal
FIG. 6 (color online). Predicted room-temperature thermal conductivity of nc-Si1−x Gex samples with x ¼ 0.5 composition and different granulometry.
conductivity of the bulk alloy is recovered. Present calculations also suggest that for hdg i < 15 nm the reduction of thermal conductivity is larger. This threshold is consistent with previous observations [22,24] and could be useful in order to tailor real samples for devices. The trend for vanishingly small grains is set by a bottom limit: in order to be stable at finite temperature, they must have a radius above the capillarity threshold [23]. According to Fig. 6, we propose a κ ∼ 1 W K−1 m−1 value as the ultimate minimum conductivity in nc-Si1−x Gex systems. In conclusion, we have calculated the thermal conductivity in nanocrystalline Si1−x Gex nanocomposites by nonequilibrium molecular dynamics simulations in a transient regime, which (although limited to effective average linear response to thermal gradients) have several appealing features, including robustness, ease of implementation, and reduced computational workload. We prove that ncSi1−x Gex systems have thermal conductivity below the corresponding bulk alloy with the same stoichiometry. The main role of affecting thermal conduction is provided by grain boundaries rather than by chemistry. Although similar interfacelike effects can be found even in short-period superlattices [25,26] or core-shell nanowires, we note that nanocomposites are admittedly superior for practical largescale commercial use since they can be obtained by lowcost bulk processing rather than by nanofabrication. Moreover, they are thermally stable systems, as required for commercial thermoelectric devices that must operate for years at high temperature. We acknowledge D. Donadio (MPI Mainz, Germany) for many helpful discussions and computational support by ISCRA-CINECA Project THEBUNA.
*
[email protected] [1] G. S. Nolas, J. Sharp, and H. Goldsmid, Thermoelectrics: Basic Principles and New Materials Developments (Springer, New York, 2001). [2] G. J. Snyder and E. S. Toberer, Nat. Mater. 7, 105 (2008).
065901-4
PRL 112, 065901 (2014)
PHYSICAL REVIEW LETTERS
[3] A. J. Minnich, M. S. Dresselhaus, Z. F. Ren, and G. Chen, Energy Environ. Sci. 2, 466 (2009). [4] M. S. Dresselhaus, G. Chen, M. Y. Tang, R. Yang, H. Lee, D. Wang, Z. Ren, and J.-P. Fleurial, Adv. Mater. 19, 1043 (2007). [5] J. H. Lienhard IV and J. H. Lienhard V, A Heat Transfer Textbook (Phlogiston Press, Cambridge, England, 2006). [6] D. A. McQuarrie, Statistical Mechanics (University Science Books, Sausalito, CA, 2000). [7] F. Müller-Plate, J. Chem. Phys. 106, 6082 (1997). [8] P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys. Rev. B 65, 144306 (2002). [9] S. G. Volz and G. Chen, Phys. Rev. B 61, 2651 (2000). [10] Y. He, I. Savić, D. Donadio, and G. Galli, Phys. Chem. Chem. Phys. 14, 16 209 (2012). [11] T. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. Dewitt, Fundamentals of Heat and Mass Transfer (John Wiley & Sons, New York, 2002). [12] E. Lampin, P. L. Palla, P.-A. Francioso, and F. Cleri, J. Appl. Phys. 114, 033525 (2013). [13] J. Garg, N. Bonini, B. Kozinsky, and N. Marzari, Phys. Rev. Lett. 106, 045901 (2011). [14] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt-Sanders, London, 1976).
week ending 14 FEBRUARY 2014
[15] J. Tersoff, Phys. Rev. Lett. 56, 632 (1986). [16] J. Tersoff, Phys. Rev. B 39, 5566 (1989). [17] P. C. Kelires and J. Tersoff, Phys. Rev. Lett. 63, 1164 (1989). [18] K. L. Whiteaker, I. K. Robinson, J. E. Van Nostrand, and D. G. Cahill, Phys. Rev. B 57, 12410 (1998). [19] S. Plimpton, J. Comput. Phys. 117, 1 (1995); see also http:// lammps.sandia.gov. [20] H. Stöhr and W. Klemm, Z. Anorg. Allg. Chem. 241, 305 (1939). [21] B. Abeles, Phys. Rev. 131, 1906 (1963). [22] X. W. Wang, H. Lee, Y. C. Lan, G. H. Zhu, G. Joshi, D. Z. Wang, J. Yang, A. J. Muto, M. Y. Tang, J. Klatsky, S. Song, M. S. Dresselhaus, G. Chen, Z. F. Ren et al., Appl. Phys. Lett. 93, 193121 (2008). [23] A. Mattoni and L. Colombo, Phys. Rev. B 78, 075408 (2008). [24] G. Joshi, H. Lee, Y. Lan, X. Wang, G. Zhu, D. Wang, R. W. Gould, D. C. Cuff, M. Y. Tang, M. S. Dresselhaus, G. Chen, and Z. Ren, Nano Lett. 8, 4670 (2008). [25] J. Garg, N. Bonini, and N. Marzari, Nano Lett. 11, 5135 (2011). [26] Z. Aksamija and I. Knezevic, Phys. Rev. B 88, 155318 (2013).
065901-5
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
Lattice Thermal Conductivity of Si1−x Gex Nanocomposites Claudio Melis and Luciano Colombo* Dipartimento di Fisica, Università di Cagliari, Cittadella Universitaria, I-09042 Monserrato (Ca), Italy (Received 24 July 2013; published 11 February 2014) We calculate the lattice thermal conductivity in model Si1−x Gex nanocomposites by molecular dynamics in a transient thermal conduction regime. Our simulations provide evidence that thermal transport depends only marginally on stoichiometry in the range 0.2 ≤ x ≤ 0.8, while it is deeply affected by the granulometry. In particular, we show that Si1−x Gex nanocomposites have lattice thermal conductivity below the corresponding bulk alloy with the same stoichiometry. The main role in affecting thermal conduction is provided by grain boundaries, which largely affect vibrational modes with a long mean-free path. DOI: 10.1103/PhysRevLett.112.065901
PACS numbers: 66.70.-f, 81.07.Bc, 44.10.+i
The physics of thermoelectric conversion is summarized by the figure of merit ZT ¼ S2 Tσ=κ characterizing the efficiency of heat conversion to electricity at temperature T (here, S is the Seebeck coefficient, while σ and κ are the electrical and thermal conductivity, respectively). Energy applications call for high ZT values or, equivalently, for minimal κ (still preserving good charge transport characteristics) [1]. In most systems of current interest, the phonon (lattice) contribution to κ is largely dominant to the electron one. This motivates the search for new materials with very low lattice thermal conductivity, a situation obtained either by tailoring new chemically complex materials with suitable phonon properties [2] or by nanostructuring more conventional materials so that internal boundaries (phase or grain boundaries, as well as interfaces) maximize phonon scattering [3,4]. In this Letter we offer a contribution to this ongoing effort by discussing, in particular, the lattice thermal conduction in nanocrystalline Si1−x Gex nanocomposites (nc-Si1−x Gex ). To this aim, we exploit a procedure to calculate κ through a MD simulation in a transient conduction (nonequilibrium) condition, which has appealing features with respect to standard techniques. By establishing a constant thermal gradient ∂T=∂z within an homogeneous and incompressible material, the corresponding steady-state thermal conduction regime is well described by the Fourier equation jz ¼ −κ∂T=∂z, where jz is the heat current [5]. By MD simulations κ is usually computed either by the Green-Kubo formalism [6] or by so-called direct methods [7,8]. In the first case, the statistical theory of transport coefficients is used which provides κ through the R current-current autocorrelation function κ ¼ ð1=kB T 2 VÞ 0þ∞ hjz ðtÞjz ð0Þidt, where kB is the Boltzmann constant, V is the system volume, and h⋅i indicates the ensemble average. Direct methods, instead, compute separately the external perturbation ∂T=∂z and the system response jz and straightforwardly provide the thermal conductivity as the response-to-perturbation ratio κ ¼ −hjz ðtÞi=∂T=∂z. Practical applications of both 0031-9007=14=112(6)=065901(5)
approaches are limited by numerical as well as conceptual restrictions. First of all, we remark that the heat current for a system of N interacting atoms is typically written as P ⃗ jðtÞ ¼ dtd Ni¼1 ϵi ðtÞ⃗ri ðtÞ, where r⃗ i and ϵi are the position and energy density of the ith particle [8,9]. Such an expression is very cumbersome to calculate except for the simplest interatomic potentials and is not directly applicable when using total-energy functionals. In addition, both hjz ðtÞjz ð0Þi and hjz ðtÞi are very much slowly convergent in time and/or size [10]. Direct methods are additionally affected by the very long simulation times needed to establish a steady-state condition, once a given temperature gradient or heat current is imposed. This often makes such methods impractical, especially for systems (as those here investigated) that, by their very structure, require large simulation cells and/or are characterized by low thermal conduction. In this work we follow a different MD approach, which simulates a transient thermal conduction regime [11,12] providing κ by directly solving the heat equation ∂T ∂ 2T ¼ κ¯ 2 ; ∂t ∂z
(1)
where κ¯ ¼ κ=ρCv is the thermal diffusivity of the system with density ρ and specific heat Cv . This method is referred to as approach-to-equilibrium MD (AEMD) [12] and it is adopted here to investigate thermal conduction in ncSi1−x Gex systems differing by stoichiometry and granulometry. We provide evidence that nc-Si1−x Gex systems have room-temperature thermal conductivity below the bulk alloy limit [13]. Furthermore, we show that κ is largely affected by granulometry, while it only marginally depends on stoichiometry over an extended range of x values. Border conditions must be imposed on Eq. (1), specifying the initial temperature profile Tðz; t ¼ 0Þ along the z direction of a simulation cell as large as Lz . Within a MD simulation the natural choice is selecting periodic
065901-1
© 2014 American Physical Society
conditions and, therefore, we assume a steplike periodic profile with Tðz; t ¼ 0Þ ¼ T 1 in the semicell 0 ≤ z < Lz =2 and Tðz; t ¼ 0Þ ¼ T 2 in the semicell Lz =2 ≤ z < Lz (with T 1 > T 2 ). It easy to prove that the general solution of Eq. (1) can be cast in the form Tðz; tÞ ¼
∞ T1 þ T2 X 2 Bn sinðαn zÞe−αn κ¯ t ; þ 2 n¼1
(2)
where αn ¼ 2πn=Lz and Bn are the coefficients of the Fourier expansion of the periodic Tðz; tÞ function. By aging the system in a microcanonical MD simulation, the initial steplike temperature profile is progressively smoothed by thermal conduction; the average temperatures hT 1 i and hT 2 i in the two semicells vary in time and the systems eventually approach a uniform temperature (i.e., equilibrium) configuration. During such a transient regime we can define the time-dependent difference in average temperatures ΔTðtÞ ¼ hT 1 ðtÞi − hT 2 ðtÞi which, through Eq. (2), is easily calculated analytically as ΔTðtÞ ¼
∞ X
2
Cn e−αn κ¯ t ;
(3)
n¼1
where Cn ¼ 8ðT 1 − T 2 Þ½cosðαn Lz =2Þ − 12 =α2n L2z . This is the key equation of AEMD: simply, the ΔTðtÞ profile must be computed on the fly during the microcanonical run and then fitted by Eq. (3) to get κ¯ . Quantum corrections to the specific heat in the corresponding thermal conductivity κ ¼ ρCv κ¯ must be duly considered below the Debye temperature of the investigated system [14]. While more details of the method are discussed elsewhere [12], we note that AEMD provides evidence that Eq. (1) is still valid at the length scales typically explored by MD simulations. Figure 1 shows the time evolution of a steplike temperature profile with ΔTð0Þ ¼ 200 K generated across a Si0.5 Ge0.5 bulk alloy sample with Lz ¼ 166.35 nm. While the smooth lines represent the 500
formal solutions of Eq. (1), the noisy lines represent the direct MD calculation of local temperature during the microcanonical run. We point out other important AEMD features, namely, as anticipated, there is no need to calculate or impose the heat current or to run long simulations to establish a steady state. The method is also very robust: no result discussed in this Letter is significantly affected by varying the initial temperature step ΔTð0Þ in the range 50–200 K or by varying the number of exponential terms in Eq. (2) in the range 10–100. In short, AEMD is characterized by a clean theoretical foundation, ease of implementation, numerical robustness, and comparatively light computational effort. By using AEMD with the Tersoff force field [15,16] (which provides good predictions on the structure of Si1−x Gex alloys [17,18], as well as reasonable results for their thermal properties [10,12]), we first calculated the room-temperature thermal conductivity of Si1−x Gex bulk alloys. In order to predict the actual values for bulk samples, for any chemical composition of the alloy, we simulated several periodically repeated simulation cells with the same lateral dimensions, 3.8 nm ≤ Lx ¼ Ly ≤ 3.9 nm (varying with the stoichiometry), and z thickness varying in the range 273 ≤ Lz ≤ 563 nm. By the standard procedure [8] of plotting 1=κ versus 1=Lz , we obtained the bulk κ value by linearly extrapolating the calculated data in the limit 1=Lz → 0. The initial steplike temperature profile was set by Nosé-Hoover thermostatting with ΔTð0Þ ¼ 200 K, while in the following microcanonical MD run the equations of motion were integrated by the velocity Verlet algorithm with a 1 fs time step. The LAMMPS code was used [19]. Our results are compared in Fig. 2 with state-of-the-art ab initio data [13] and available experiments [20–22]. The resulting picture represents the reliability of the present method: κ is basically independent of the stoichiometry in the range 0.2 ≤ x ≤ 0.8, but its value sharply increases toward the c-Si and c-Ge values as x → 0 or 1, respectively. In particular, we obtained 216 18 and 65 6 W K−1 m−1 250
t=0 ps t=20 ps t=50 ps t=100 ps
450
This work Ref. [13] Ref. [20] Ref. [21] Ref. [22]
200
400
-1 -1 κ [W K m ]
Temperature [K]
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
PRL 112, 065901 (2014)
350 300 250
150 100 50
200 150
0 0
20
40
60
80 100 120 140 160
0
z [nm]
0.2
0.4
0.6
0.8
1
% Germanium Content
FIG. 1 (color online). Evolution of the temperature profile across a Si0.5 Ge0.5 bulk alloy sample with Lz ¼ 166.35 nm as provided by solving analytically Eq. (1) (smooth lines) and by direct MD calculation of local temperature (noisy color lines).
FIG. 2 (color online). Room-temperature thermal conductivity of bulk Si1−x Gex alloys predicted by AEMD simulations (red squares). Experimental data (color triangles and diamonds) and ab initio calculations (black dots) are shown for comparison.
065901-2
PRL 112, 065901 (2014)
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
for c-Si and c-Ge to compare with the corresponding range of experimental values 130–240 and 50–60 W K−1 m−1 , respectively (see Ref. [10] and references therein). While the agreement for Ge is very good, we remark that our result for Si is better than most of the previously reported values obtained with other methods, by using the same Tersoff force field. Only the calculation by He et al. [10] obtained a slightly better κ for c-Si, using, however, a huge simulation cell and extremely long simulation times. Generating atomistic models of nc-Si1−x Gex is a challenging and very computer-intensive problem, since we need to correctly reproduce both local-scale and global structural features (respectively, grain boundary atomic architecture and grain size or shape distribution) for many different stoichiometries. With this in mind, we followed a multistep procedure repeatedly applied to simulation cells that, by construction, initially contained just Si atoms and had lateral dimensions Lx ¼ 2.715 nm and Ly ¼ 27.15 nm. The z thickness was instead varied in the range 54.3 ≤ Lz ≤ 135.7 nm for the reason discussed above. Each system was at first fully amorphized by quenching from the melt. Then, N g sites were selected at random in the yz plane and around each of them a cylindrical void (passing across the full Lx dimension) was created by removing atoms. Voids were distributed avoiding overlap, and their initial radius was randomly set above the capillarity threshold [23] (∼0.5 nm). Next, each void was filled by a crystalline cylinder randomly rotated in the yz plane. The resulting mixed structures (an amorphous matrix spotted by crystalline seeds) were annealed at constant T ¼ 1200 K for very long time, varying from 1.5 to 4.0 ns for increasing granularity. During this step seed grains underwent growth until nanocrystalline silicon systems (hereafter referred to as templates) were eventually obtained. Three different templates are shown in Fig. 3: the resulting granularity hdg i reported in the figure for each sample
was effectively controlled by the initial choice of N g . At this stage nc-Si1−x Gex samples were obtained by randomly replacing Si atoms with Ge atoms in the templates, followed by further careful relaxation. This alchemy was flanked by a self-affine rescaling of the atomic positions, so that the actual lattice spacing of any nc-Si1−x Gex sample was set as aðSi1−x Gex Þ ¼ ð1 − xÞaðSiÞ þ xaðGeÞ, where aðSiÞ and aðGeÞ are the Si and Ge lattice constant provided by the Tersoff potential, respectively. This somewhat complex procedure is not only effective in producing trustworthy atomistic models (compare Fig. 3 with TEM pictures reported in Refs. [22,24]), but it is also very efficient: as a matter of fact, the most computerintensive step (i.e., the growth of crystalline grains against amorphous regions) must be followed just once, while Si → Ge replacements and following relaxation are comparatively light and, therefore, can be effortlessly repeated for any stoichiometry we are interested in. In total we generated a set of five different templates, with 3 ≤ hdg i ≤ 25 nm, corresponding to a number of atoms varying in the range ð2–5Þ × 105 (according to the actual Lz value of the simulation cell), indeed a large enough number to warrant that the random chemical arrangement of Si and Ge atoms is representative of real samples at any investigated stoichiometry. More importantly, however, the explored range of average grain size perfectly matches real nc-Si1−x Gex samples [22,24]. The values of thermal conductivity (for any given stoichiometry and granularity) discussed below were obtained as the 1=Lz → 0 extrapolation of a set of κ values corresponding to the five different systems previously generated. Furthermore, the thermal conductivity of each single system was calculated as the configurational average over two different samples with the same Lz , x, and hdg i. We first address the key question of whether thermal conductivity is larger in bulk alloys or in nanocomposites. In Fig. 4 the remarkable result of the present investigation is reported for a nc-Si1−x Gex sample of granularity hdg i ¼ 10.5 nm. Here we provide evidence that κ is smaller by a ∼50% in nanocomposites, consistent with experimental data which report a very similar reduction in thermal 6
Bulk alloy Nanocomposite
κ [W K -1 m-1]
5 4 3 2 1
FIG. 3 (color online). Left: Snapshots of the final configuration obtained after the growth process described in the text for three template samples differing for average grain size hdg i. Right: Enlargement of the marked region, showing the atomistic structure of grain boundaries.
0.2
0.3
0.4
0.5
0.6
0.7
0.8
% Germanium Content
FIG. 4 (color online). Predicted room-temperature thermal conductivity of nc-Si1−x Gex samples with hdg i ¼ 10.5 nm and different stoichiometry.
065901-3
90
4 Bulk alloy Nanocomposite
80
3
70 60 50 40 30
2.5 2 1.5 1
20
0.5
10 0 0.0001 0.001
Nanocomposite Bulk alloy
3.5
-1 κ [W K m-1]
Thermal conductivity accumulation [%]
100
week ending 14 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
PRL 112, 065901 (2014)
0.01
0.1
1
0
10
0
5
10
15
20
25
30
Average grain size [nm]
Lz [µm]
FIG. 5 (color online). Thermal conductivity accumulation function for a nc-Si0.5 Ge0.5 sample (black dots) and a Si0.5 Ge0.5 bulk alloy (red squares). Solid lines show the corresponding accumulation extracted from the linear 1=κ versus 1=Lz extrapolation.
conductivity of both p- and n-type samples [22,24]. While experiments were collected for the sole x ¼ 0.2 stoichiometry, Fig. 4 shows that this feature is indeed valid for any chemical composition. This is an important evidence supporting the use of nc-Si1−x Gex systems as thermoelectric materials of possibly large technological impact. The reduced value of κ is clearly due to the enhancement of scattering phenomena provided by the network of grain boundaries (see Fig. 3). This statement is consistent with recent results on short-period Si=Ge superlattices [25,26], where the major role played by interface scattering has been clearly defined. We further substantiate this concept by the calculation of the thermal conductivity accumulation function, here defined as the ratio κðLz Þ=κ ∞ between its value computed for a simulation cell with length Lz and its extrapolated value for a bulklike (infinite) sample. κðLz Þ=κ ∞ gives the contribution to the thermal conductivity provided by vibrational modes with mean-free path λ up to Lz . In Fig. 5 it is shown that for a nc-Si0.5 Ge0.5 sample, as much as ∼80% of the thermal conductivity is provided by vibrational modes with λ ≤ 100 nm. On the other hand, in the bulk alloy the same contribution is provided by vibrations with much longer λ, up to 1 μm, consistent with the estimation by Garg et al. [13]. Another feature emerging from Fig. 4 is that, similarly to bulk alloys, the actual chemical composition seems to play a minor role: as a matter of fact, κ is basically independent of stoichiometry, although there is an indication that a small Ge content should be preferred in order to fully minimize κ. We can then draw the conclusion that alloying is not the ultimate solution for minimizing thermal conduction in Si1−x Gex systems; rather, this role is played by grain boundaries. In order to better characterize the role of the nanostructure, we have extended the above calculations to the full set of nc-Si1−x Gex samples with x ¼ 0.5. Figure 6 shows that κ increases monotonically with granulometry, and, for a large enough average grain size, the thermal
FIG. 6 (color online). Predicted room-temperature thermal conductivity of nc-Si1−x Gex samples with x ¼ 0.5 composition and different granulometry.
conductivity of the bulk alloy is recovered. Present calculations also suggest that for hdg i < 15 nm the reduction of thermal conductivity is larger. This threshold is consistent with previous observations [22,24] and could be useful in order to tailor real samples for devices. The trend for vanishingly small grains is set by a bottom limit: in order to be stable at finite temperature, they must have a radius above the capillarity threshold [23]. According to Fig. 6, we propose a κ ∼ 1 W K−1 m−1 value as the ultimate minimum conductivity in nc-Si1−x Gex systems. In conclusion, we have calculated the thermal conductivity in nanocrystalline Si1−x Gex nanocomposites by nonequilibrium molecular dynamics simulations in a transient regime, which (although limited to effective average linear response to thermal gradients) have several appealing features, including robustness, ease of implementation, and reduced computational workload. We prove that ncSi1−x Gex systems have thermal conductivity below the corresponding bulk alloy with the same stoichiometry. The main role of affecting thermal conduction is provided by grain boundaries rather than by chemistry. Although similar interfacelike effects can be found even in short-period superlattices [25,26] or core-shell nanowires, we note that nanocomposites are admittedly superior for practical largescale commercial use since they can be obtained by lowcost bulk processing rather than by nanofabrication. Moreover, they are thermally stable systems, as required for commercial thermoelectric devices that must operate for years at high temperature. We acknowledge D. Donadio (MPI Mainz, Germany) for many helpful discussions and computational support by ISCRA-CINECA Project THEBUNA.
*
[email protected] [1] G. S. Nolas, J. Sharp, and H. Goldsmid, Thermoelectrics: Basic Principles and New Materials Developments (Springer, New York, 2001). [2] G. J. Snyder and E. S. Toberer, Nat. Mater. 7, 105 (2008).
065901-4
PRL 112, 065901 (2014)
PHYSICAL REVIEW LETTERS
[3] A. J. Minnich, M. S. Dresselhaus, Z. F. Ren, and G. Chen, Energy Environ. Sci. 2, 466 (2009). [4] M. S. Dresselhaus, G. Chen, M. Y. Tang, R. Yang, H. Lee, D. Wang, Z. Ren, and J.-P. Fleurial, Adv. Mater. 19, 1043 (2007). [5] J. H. Lienhard IV and J. H. Lienhard V, A Heat Transfer Textbook (Phlogiston Press, Cambridge, England, 2006). [6] D. A. McQuarrie, Statistical Mechanics (University Science Books, Sausalito, CA, 2000). [7] F. Müller-Plate, J. Chem. Phys. 106, 6082 (1997). [8] P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys. Rev. B 65, 144306 (2002). [9] S. G. Volz and G. Chen, Phys. Rev. B 61, 2651 (2000). [10] Y. He, I. Savić, D. Donadio, and G. Galli, Phys. Chem. Chem. Phys. 14, 16 209 (2012). [11] T. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. Dewitt, Fundamentals of Heat and Mass Transfer (John Wiley & Sons, New York, 2002). [12] E. Lampin, P. L. Palla, P.-A. Francioso, and F. Cleri, J. Appl. Phys. 114, 033525 (2013). [13] J. Garg, N. Bonini, B. Kozinsky, and N. Marzari, Phys. Rev. Lett. 106, 045901 (2011). [14] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt-Sanders, London, 1976).
week ending 14 FEBRUARY 2014
[15] J. Tersoff, Phys. Rev. Lett. 56, 632 (1986). [16] J. Tersoff, Phys. Rev. B 39, 5566 (1989). [17] P. C. Kelires and J. Tersoff, Phys. Rev. Lett. 63, 1164 (1989). [18] K. L. Whiteaker, I. K. Robinson, J. E. Van Nostrand, and D. G. Cahill, Phys. Rev. B 57, 12410 (1998). [19] S. Plimpton, J. Comput. Phys. 117, 1 (1995); see also http:// lammps.sandia.gov. [20] H. Stöhr and W. Klemm, Z. Anorg. Allg. Chem. 241, 305 (1939). [21] B. Abeles, Phys. Rev. 131, 1906 (1963). [22] X. W. Wang, H. Lee, Y. C. Lan, G. H. Zhu, G. Joshi, D. Z. Wang, J. Yang, A. J. Muto, M. Y. Tang, J. Klatsky, S. Song, M. S. Dresselhaus, G. Chen, Z. F. Ren et al., Appl. Phys. Lett. 93, 193121 (2008). [23] A. Mattoni and L. Colombo, Phys. Rev. B 78, 075408 (2008). [24] G. Joshi, H. Lee, Y. Lan, X. Wang, G. Zhu, D. Wang, R. W. Gould, D. C. Cuff, M. Y. Tang, M. S. Dresselhaus, G. Chen, and Z. Ren, Nano Lett. 8, 4670 (2008). [25] J. Garg, N. Bonini, and N. Marzari, Nano Lett. 11, 5135 (2011). [26] Z. Aksamija and I. Knezevic, Phys. Rev. B 88, 155318 (2013).
065901-5
