Sub-amorphous thermal conductivity in amorphous heterogeneous nanocomposites

Abstract

Pure amorphous solids are traditionally considered to set the lower bound of thermal conductivity due to their disordered atomic structure that impedes vibrational energy transport. However, the lower limits for thermal conductivity in heterogeneous amorphous solids and the physical mechanisms underlying these limits remain unclear. Here, we use equilibrium molecular dynamics to show that an amorphous SiGe nanocomposite can possess thermal conductivity substantially lower than those of the amorphous Si and Ge constituents. Normal mode analysis indicates that the presence of the Ge inclusion localizes vibrational modes with frequency above the Ge cutoff in the Si host, drastically reducing their ability to transport heat. This observation suggests a general route to achieve exceptionally low thermal conductivity in fully dense solids by restricting the vibrational density of states available for transport in heterogeneous amorphous nanocomposites.

---

1 Introduction

Low thermal conductivity materials are desired for a wide range of applications ranging from thermoelectric power generators^1–6 to thermopile detectors.⁷ Traditionally, amorphous materials are considered to set the lower limit of thermal conductivity due to the disordered atomic structure that impedes the formation of propagating vibrations.^8–10 While in crystals heat is carried by propagating lattice waves, or phonons, in amorphous solids the lack of a periodic atomic structure results in very different mechanisms for vibrational energy transport.

Allen and Feldman introduced categories of vibrational modes in amorphous solids known as propagons, diffusons, and locons.^11,12 Propagons are propagating and delocalized phonon-like plane waves that typically possess long wavelengths compared to the interatomic spacing. Diffusons are modes that scatter over a distance less than their wavelength and thus transport heat as a random-walk. Locons are non-propagating and localized modes that are unable to transport heat in harmonic solids.^12,13

This classification has been widely used to interpret experiments and calculations of transport in amorphous materials, particularly for pure a-Si. For instance, numerical works using equilibrium molecular dynamics (EMD) and lattice dynamics (LD) have attempted to determine the fraction of heat carried by each type of vibration. In their original work, Allen et al. reported that ∼20% of thermal conductivity of a-Si is from propagons (≲3 THz) whereas the rest are from diffusons (3–17 THz) and none is from locons (≳17 THz).¹⁴ He et al. reported that although only 3% of the mode population is propagons, they transport up to 50% of the heat due to their long propagation distances.¹⁵ Calculations by Larkin and McGaughey indicate that propagons have a lifetime scaling of ω⁻² which suggests that these modes are plane-wave-like and are propagating.¹³ Wei and Henry have reported that frequency modes above ∼17 THz are highly localized and do not contribute to thermal conductivity using Green-Kubo modal analysis for a-Si.¹⁶

Experimental works have qualitatively confirmed some of these predictions.^17–19 Sultan et al. reported that modification of the surface of an amorphous SiN membrane changes the thermal conductance of the membrane, indicating the importance of propagons for heat conduction.²⁰ They estimated that propagons are responsible for ∼40–50% of thermal conductivity in amorphous SiN using kinetic theory. Braun et al. reported that diffusons are the dominant heat carriers for films of thickness less than 100 nm, while the propagon contribution is present in thicker films.¹⁹

Although pure amorphous solids are typically assumed to achieve the lower limit of thermal conductivity, some works have examined how this limit may be broken. In semi-crystalline solids, it is well known that composites can possess exceptionally low thermal conductivity due to thermal boundary resistance.²¹ This effect has been exploited by Chiritescu et al.²² to achieve ultralow thermal conductivity in disordered WSe₂ nanolaminates below the minimum thermal conductivity predicted by the Cahill–Pohl model,²³ although a recent theory work suggests that the experiments agree with this model if anisotropy is taken into account.²⁴ Wingert et al. reported that crystalline silicon nanotubes with shell thicknesses as thin as 5 nm have a low thermal conductivity of 1.1 W m⁻¹ K⁻¹, lower than that of the amorphous counterpart via a phonon softening effect.²⁵ Dechaumphai et al. experimentally observed an ultralow thermal conductivity of 0.33 ± 0.04 W m⁻¹ K⁻¹ at room temperature in amorphous multilayers made of Au and Si.²⁶ Computationally, Norouzzadeh et al. used MD to study the thermal conductivity of an a-SiGe alloy with different Ge content and observed thermal conductivity values below those of the constituent materials.²⁷ Giri et al. used NEMD to examine the role of the interface of amorphous SiGe superlattices and amorphous Si/heavy-Si superlattices, concluding that increasing mass-mismatch in amorphous superlattices results in higher Kapitza resistances, leading to low thermal conductivity.²⁸

Although these works have suggested that thermal conductivities of heterogeneous amorphous solids below those of the pure constituents are achievable, key questions remain. Some of these works have interpreted their results with a phonon gas model, which is of questionable validity for diffusons and locons, and others have used the concept of thermal boundary resistance to explain their observations. In particular, the latter approach implicitly assumes that vibrational modes of the two solids composing the interface are well defined. However, if the inclusion in the nanocomposite is sufficiently small, the vibrational modes of the composite may not coincide with the vibrations of the pure materials. In this case, the nature of the vibrations in the composite solids and hence the lower limits of thermal conductivity in heterogeneous amorphous solids remain unclear.

Here, we examine heat transport in amorphous SiGe nanocomposites consisting of a Ge inclusion in a Si host matrix. We find that these structures can possess thermal conductivities that are significantly smaller than those of the constituent materials, with the minimum thermal conductivity reaching as low as 32% of that of the amorphous Si host. Lattice dynamics analysis demonstrates that the presence of the Ge cluster drastically enhances localization of vibrational modes with frequency above the Ge cutoff in the Si host, leading to a remarkable decrease in thermal conductivity. These results demonstrate a mechanism for achieving remarkably low thermal conductivity in fully dense amorphous materials that may be useful for solid-state thermal insulation and highly sensitive thermopile detectors.

2 Methods

We calculated the thermal conductivity of amorphous Si and amorphous SiGe nanocomposites using equilibrium MD with the Stillinger–Weber (SW) interatomic potential.²⁹ The two types of structures studied are shown in Fig. 1. The atomic configuration consisting of 4096 atoms was provided by N. Mousseau and was generated from the modified Wooten–Winer–Weaire (WWW) algorithm.³⁰ For na-SiGe structures, a cubic domain in the middle of the structure with side length a was replaced with heavier germanium atoms with appropriate coefficient changes in SW potential. SW potential coefficients for silicon and germanium interactions are described in ref. 29, 31 and 32 The side length, a, was chosen to be 10, 20, 25, 30, 35, 40, and 46.4 Å. These lengths represent 1, 10, 20, 35, 55, 82, and 100% Ge fraction, respectively. Periodic boundary conditions were imposed for all the structures. The MD simulations were performed with Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)³³ with a time step of 0.5 fs. The simulation procedure began with an anneal at 1000 K for 20 ns using the NPT ensemble to reduce metastabilities.^13,15 We observed a decrease and plateau of the potential energy during the annealing process for each structure indicating a reduction of metastability.

Fig. 1 4096-atom configurations of (a) amorphous silicon and (b) nanostructured amorphous silicon germanium. Blue atoms represent silicon and orange atoms represent germanium. The germanium cubic side length, a, varies from 10 Å to the side length of the entire domain, L = 46.4 Å.

Subsequently, the domain was quenched at a rate of 10 K ps⁻¹ to 300 K and equilibrated in an NPT ensemble at 300 K for 20 ns to relax the structure to equilibrium pressure. Because volume and pressure fluctuate in MD simulations, we computed the average atom positions over the last 100 ps to ensure the domain was not under strain. The resulting mean pressure was on the order of 0.1 bar. This domain was then thermostatted in an NVT ensemble for 10 ns using a Nose–Hoover thermostat. After an additional NVE equilibration for 50 ps, the heat fluxes were computed for 1.6 ns in NVE ensemble.

We computed the thermal conductivity of the various structures using the Green–Kubo (GK) formalism, which relates the thermal conductivity to the heat current autocorrelation function by

where k_B is the Boltzmann constant, T is the temperature, V is the system volume, t is time, and J is the heat flux. The angular brackets denote an ensemble average. The thermal conductivity calculations reported in this study are based on the average of the integrals of the heat current autocorrelation functions (HCACF) from 10 simulations.

Fig. 2(a) shows the HCACF normalized by 〈J(0)·J(0)〉 for a-Si. The autocorrelation function converges quickly to 0 in less than 0.5 ps for a-Si. The HCACF convergence times of na-SiGe are on the order of 10 ps. The resulting thermal conductivity of a-Si obtained from the integral of the autocorrelation function versus integration time is depicted in Fig. 2(b). The thermal conductivity of a-Si is determined by taking the average between 5 and 20 ps. The thermal conductivity of a-Si with respect to temperature for 4096 atoms with SW potential is plotted in Fig. 2(c) and compared with works by Larkin and McGaughey¹³ and Lv and Henry.¹⁶ At 300 K, thermal conductivity from this work is 1.55 ± 0.20 W m⁻¹ K⁻¹ which is in agreement with these works. Consistent with ref. 16, weak temperature dependence of thermal conductivity is observed. Direct comparison to experimental results is difficult as thermal conductivity of a-Si varies significantly by the fabrication process, hydrogenation, heat treatment, and defects, but experimental thermal conductivity typically ranges from 1 to 6 W m⁻¹ K⁻¹ at room temperature.^15,17,18,34

Fig. 2 (a) Normalized heat current autocorrelation function versus time. (b) Thermal conductivity integral calculated by eqn (1)versus time. The thermal conductivity of a-Si is determined by taking the average between 5 and 20 ps. (c) Thermal conductivity versus temperature (blue circles), comparison with the works by Larkin and McGaughey (black squares),¹³ and Lv and Henry (red crosses)¹⁶ utilizing 4096 atoms, SW potential, and GK formalism at temperatures from 300 K to 1000 K. No temperature dependence is observed.

3 Results

We now examine the thermal conductivity of na-SiGe versus Ge content, shown in Fig. 3. Pure amorphous Si and Ge have thermal conductivities of 1.55 ± 0.20 W m⁻¹ K⁻¹ and 0.99 ± 0.21 W m⁻¹ K⁻¹, respectively. Interestingly, we observe thermal conductivities substantially smaller than either of these values for na-SiGe composites with Ge content ranging from 35% to 82%, with the minimum thermal conductivity of 0.50 ± 0.17 W m⁻¹ K⁻¹ achieved with 55% of Ge content. This value is less than a third of the original a-Si thermal conductivity and half that of a-Ge. Interestingly, it is even lower than the thermal conductivity of an amorphous SiGe alloy with the same Ge content, which has a thermal conductivity of 0.78 ± 0.16 W m⁻¹ K⁻¹. The percentage decrease of thermal conductivity in na-SiGe is nearly twice that in a-Si/a-Ge superlattices by an NEMD study by Giri et al. utilizing SW potential despite similar geometry.²⁸

Fig. 3 Thermal conductivity of na-SiGe versus Ge content. The minimum thermal conductivity of 0.50 ± 0.17 W m⁻¹ K⁻¹ is observed with 55% Ge content.

To understand the mechanism behind the reduction in thermal conductivity, we first examine the vibrational density of states (vDOS) of the a-Si and a-Ge constituents in na-Si_0.45Ge_0.55 shown in Fig. 4(a). The vDOS is computed from

where N_atom is the number of atoms, T is the temperature, m_n is the mass of atom n, and v_n(t) is the velocity of atom n at time t.³⁵ The vDOS of a-Si and a-Ge is similar to that of c-Si and c-Ge with distinct peaks at certain frequencies.²⁸ Due to absence of strong anharmonicity, only weak vibrational interaction of Si and Ge atoms is expected for frequencies greater than the frequency cutoff of a-Ge of 10 THz. In other words, we expect the vibrational modes with frequencies exceeding 10 THz to be confined to a-Si.

Fig. 4 (a) The vibrational density of states of Si (blue dashed line) and Ge (red dashed line) constituents in na-Si_0.45Ge_0.55 along with the total density of states (black line). Inverse participation ratio (IPR) for (b) a-Si, (c) na-Si_0.90Ge_0.10, and (d) na-Si_0.45Ge_0.55. (e) Zoomed-in view of IPR of na-Si_0.45Ge_0.55 for frequencies from 5 to 15 THz where modes above the bold line are defined as locons. Vibrational modes start to be localized at 9 THz and are completely localized above 10 THz.

We confirm this hypothesis by first calculating the inverse participation ratio (IPR), which is a measure of how many atoms participate in the motion of a particular eigenmode. The IPR is given by

where e_iα,n is the eigenvector component for atom i in α direction for the mode n.³⁶ The eigenvectors for each mode and atom are calculated by harmonic lattice dynamics in GULP³⁷ with relaxed structures from MD at 300 K. The IPR is defined so that it equals 1/N_atom if all atoms are participating, or 1 if the vibration is completely localized to one atom. Defining a specific IPR value that uniquely distinguishes locons is not possible, but vibrational modes with participation ratio less than 0.2 (corresponding to IPR greater than 0.0012 here) have been defined previously as localized modes.^38,39 We therefore define locons according to this convention.

Fig. 4(b)–(e) show the IPR for a-Si, na-Si_0.90Ge_0.10, na-Si_0.45Ge_0.55, and a zoomed-in view of the IPR of na-Si_0.45Ge_0.55 from 5 to 15 THz. The IPR for a-Si, Fig. 4(b), shows that locons are observed primarily over around 17 THz, consistent with prior works.^12,16 As Ge atoms are introduced in the nanocomposite in na-Si_0.90Ge_0.10, we observe locons in the medium-frequency region around 10 THz. For na-Si_0.45Ge_0.55, all the vibrational modes above around 10 THz are localized. The corresponding locon mode fractions are 7%, 9%, and 31% for a-Si, na-Si_0.90Ge_0.10, and na-Si_0.45Ge_0.55, respectively. In other words, na-Si_0.45Ge_0.55 has the lowest thermal conductivity and also more than 4 times the number of locons than a-Si, suggesting localized modes in Si are associated with the low thermal conductivity of the nanocomposite. We also note that vibrational modes with higher IPR than 0.0012 are present at low frequencies. We have verified that these modes are due to the finite size of the computational domain and disappear as the size of the system increases.

We next confirm that these localized modes reside in silicon by calculating the local vibrational density of states, defined as⁴⁰

where the sum is over Cartesian directions α and vibrational modes n for atom i. Furthermore, the spatial distribution of energy can be described as³⁹where n_BE is the occupation number given by the Bose–Einstein distribution. We identify where the vibrational modes are localized by performing the sum only for vibrational modes that correspond to locons as identified by the IPR.

The spatial energy distribution is shown in Fig. 5(a)–(c) for a-Si, na-Si_0.90Ge_0.10, and na-Si_0.45Ge_0.55, respectively. The distribution has been normalized by the maximum energy of an atom in the domain. We plot cross section x–y plane in the middle of z axis for clear visualization. It is apparent that for a-Si the spatial distribution of locons is randomly distributed. As Ge content is increased, however, we observe that locons are located in Si atoms. This result confirms that vibrational modes over around 10 THz are increasingly localized as Ge content grows and that these locons are indeed localized in a-Si atoms.

Fig. 5 Normalized spatial energy distribution of the cross section xy plane in the middle of z axis for (a) a-Si, (b) na-Si_0.90Ge_0.10, and (c) na-Si_0.45Ge_0.55. Individual circles in the figure represent atoms and dashed lines represent the boundaries between Si and Ge atoms. Color indicates the degree of localization at an atom with red indicating high localization. (d) Spectral thermal diffusivities of a-Si and na-Si_0.45Ge_0.55versus mode frequency. Thermal diffusivities decrease significantly for vibrational modes with frequencies higher than 10 THz in na-Si_0.45Ge_0.55 compared to those in a-Si.

The drastic increase in locon population in na-SiGe suggests that the origin of the low thermal conductivity in na-SiGe is due to conversion of delocalized modes in a-Si to locons. To verify this hypothesis, we calculate the thermal diffusivities using the harmonic heat flux operator. The thermal conductivity of a solid is given by

where V is the volume of structure, C (ω_i) is the specific heat, D_th (ω_i) is the thermal diffusivity of frequency ω_i, and the summation is over modes. For diffusons under harmonic Allen-Feldman (AF) theory, the thermal diffusivity is calculated bywhere S_ij is the heat current operator in the harmonic approximation.⁴¹ Feldman et al. demonstrated that diffusivity calculations based on Peierls–Boltzmann theory (phonon gas model) for low frequency propagon modes coincide reasonably well with D_AF in the low frequency range.⁴⁰ Therefore, we calculate D_AF for all the vibrational modes for a-Si and na-Si_0.45Ge_0.55 as shown in Fig. 5(d). We observe that, for vibrational modes defined as locons by IPR (≳17 THz for a-Si and ≳10 THz for na-Si_0.45Ge_0.55), the thermal diffusivities decrease significantly. For vibrational modes with frequencies between 10 THz to 17 THz, we observe an order of magnitude decrease in diffusivity from a-Si to na-Si_0.45Ge_0.55, contributing to the decrease in thermal conductivity. For low frequency propagating modes ≲2 THz, no apparent changes in diffusivity occur among different structures, and we are unable to conclude how propagons with frequencies less than 1 THz are affected in the nanocomposite due to limitations in the size of the domain.

The results suggest a simple explanation for the low thermal conductivity of the nanocomposite. In a-Si, nearly the full vibrational spectrum contributes to heat conduction as indicated by the calculated thermal diffusivities and associated small locon population. In the nanocomposite, diffusons with frequencies above the soft Ge cutoff frequency become localized, impeding their ability to transport heat. In effect, the soft inclusion restricts the vibrational spectrum available to conduct heat because many Si vibrational modes are not supported in the inclusion.

Another interesting consideration is why the thermal conductivity of the composite is less than the intrinsic thermal conductivity of the amorphous Ge. The explanation can again be identified from the locon population and the density of vibrational states. The locon population of a-Ge is found to be 8% compared to 31% of na-Si_0.45Ge_0.55, which means there are more vibrational modes that are able to transport heat than in na-Si_0.45Ge_0.55. Although the a-Ge has a lower cutoff frequency than a-Si, its density of states is the same as that as a-Si because the atomic number densities are identical. However, in the nanocomposite, only a fraction of the modes in Si with frequencies below that of the Ge cutoff are able to conduct heat; therefore, the composite contains fewer states with non-negligible thermal diffusivities than a-Ge. As a result, the thermal conductivity of the composite may be lower than those of both the stiff host and softer inclusion.

Many prior works have interpreted thermal conductivity reductions in amorphous or disordered heterogeneous solids using the concept of thermal boundary resistance between the adjacent layers.^4,26,28 However, this interpretation relies on the vibrational mode properties of individual constituents separately. Our analysis shows that the vibration mode characters change drastically from a-Si to na-SiGe, suggesting that thermal boundary resistance is not a well-defined concept in the amorphous nanocomposites studied here as the vibrational modes of the constituent materials cannot be separated. Instead, it is the change in character of the overall vibrational modes of the composite that leads to the low thermal conductivity.

4 Conclusions

In summary, we have studied thermal transport in amorphous heterogeneous nanocomposites using molecular dynamics and lattice dynamics. We find that the thermal conductivity of na-Si_0.45Ge_0.55 is substantially lower than that of both constituent materials due to the localization of vibrational modes in the stiff host a-Si with frequencies exceeding the cutoff of the soft inclusion. This observation suggests a general route to achieve exceptionally low thermal conductivity in fully dense amorphous solids by restricting the vibrational density of states for transport in heterogeneous nanocomposites.

Acknowledgements

This work was supported by the Samsung Scholarship, NSF CAREER Award CBET 1254213, and Boeing under Boeing-Caltech Strategic Research and Development Relationship Agreement. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. The authors thank Normand Mousseau for providing the atomic positions from the WWW algorithm, and Andrew Robbins, Benoit Latour, Wei Lv, and Asegun Henry for useful discussions.

References

1. B. Poudel, Q. Hao, Y. Ma, Y. Lan, A. Minnich, B. Yu, X. Yan, D. Wang, A. Muto, D. Vashaee, X. Chen, J. Liu, M. S. Dresselhaus, G. Chen and Z. Ren, Science, 2008, 320, 634–638.
2. K. Biswas, J. He, I. D. Blum, C.-I. Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid and M. G. Kanatzidis, Nature, 2012, 489, 414–418.
3. M. Zebarjadi, K. Esfarjani, M. S. Dresselhaus, Z. F. Ren and G. Chen, Energy Environ. Sci., 2012, 5, 5147–5162.
4. D. G. Cahill, P. V. Braun, G. Chen, D. R. Clarke, S. Fan, K. E. Goodson, P. Keblinski, W. P. King, G. D. Mahan, A. Majumdar, H. J. Maris, S. R. Phillpot, E. Pop and L. Shi, Appl. Phys. Rev., 2014, 1, 011305.
5. L.-D. Zhao, V. P. Dravid and M. G. Kanatzidis, Energy Environ. Sci., 2014, 7, 251–268.
6. A. J. Minnich, J. Phys.: Condens. Matter, 2015, 27, 053202.
7. M. C. Foote and E. W. Jones, Aerospace/Defense Sensing and Controls, 1998, pp. 192–197.
8. D. G. Cahill and R. O. Pohl, Phys. Rev. B: Condens. Matter Mater. Phys., 1987, 35, 4067.
9. A. McGaughey and M. Kaviany, Advances in Heat Transfer, Elsevier, 2006, vol. 39, pp. 169–255.
10. H. Mizuno, S. Mossa and J.-L. Barrat, Sci. Rep., 2015, 5, 14116.
11. P. B. Allen and J. L. Feldman, Phys. Rev. Lett., 1989, 62, 645.
12. P. B. Allen, J. L. Feldman, J. Fabian and F. Wooten, Philos. Mag. B, 1999, 79, 1715–1731.
13. J. M. Larkin and A. J. H. McGaughey, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 144303.
14. J. L. Feldman, M. D. Kluge, P. B. Allen and F. Wooten, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 12589.
15. Y. He, D. Donadio and G. Galli, Appl. Phys. Lett., 2011, 98, 144101.
16. W. Lv and A. Henry, New J. Phys., 2016, 18, 013028.
17. D. G. Cahill, M. Katiyar and J. R. Abelson, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 6077.
18. B. L. Zink, R. Pietri and F. Hellman, Phys. Rev. Lett., 2006, 96, 055902.
19. J. L. Braun, C. H. Baker, A. Giri, M. Elahi, K. Artyushkova, T. E. Beechem, P. M. Norris, Z. C. Leseman, J. T. Gaskins and P. E. Hopkins, Phys. Rev. B, 2016, 93, 140201.
20. R. Sultan, A. D. Avery, J. M. Underwood, S. J. Mason, D. Bassett and B. L. Zink, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 214305.
21. A. Minnich and G. Chen, Appl. Phys. Lett., 2007, 91, 073105.
22. C. Chiritescu, D. G. Cahill, N. Nguyen, D. Johnson, A. Bodapati, P. Keblinski and P. Zschack, Science, 2007, 315, 351–353.
23. D. G. Cahill, S. K. Watson and R. O. Pohl, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6131–6140.
24. Z. Chen and C. Dames, Appl. Phys. Lett., 2015, 107, 193104.
25. M. C. Wingert, S. Kwon, M. Hu, D. Poulikakos, J. Xiang and R. Chen, Nano Lett., 2015, 15, 2605–2611.
26. E. Dechaumphai, D. Lu, J. J. Kan, J. Moon, E. E. Fullerton, Z. Liu and R. Chen, Nano Lett., 2014, 14, 2448–2455.
27. P. Norouzzadeh, A. Nozariasbmarz, J. S. Krasinski and D. Vashaee, J. Appl. Phys., 2015, 117, 214303.
28. A. Giri, P. E. Hopkins, J. G. Wessel and J. C. Duda, J. Appl. Phys., 2015, 118, 165303.
29. F. H. Stillinger and T. A. Weber, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 5262.
30. G. T. Barkema and N. Mousseau, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 4985.
31. K. Ding and H. C. Andersen, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 6987.
32. M. Laradji, D. P. Landau and B. Dünweg, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 4894.
33. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
34. X. Liu, J. L. Feldman, D. G. Cahill, R. S. Crandall, N. Bernstein, D. M. Photiadis, M. J. Mehl and D. A. Papaconstantopoulos, Phys. Rev. Lett., 2009, 102, 035901.
35. J. M. Dickey and A. Paskin, Phys. Rev., 1969, 188, 1407.
36. R. J. Bell and P. Dean, Discuss. Faraday Soc., 1970, 50, 55–61.
37. J. D. Gale, J. Chem. Soc., Faraday Trans., 1997, 93, 629–637.
38. L. Yang, N. Yang and B. Li, Nano Lett., 2014, 14, 1734–1738.
39. J. Chen, G. Zhang and B. Li, Nano Lett., 2010, 10, 3978–3983.
40. J. Feldman and N. Bernstein, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 235214.
41. P. B. Allen and J. L. Feldman, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 12581.

---