View PDF VersionPrevious ArticleNext Article
 
DOI: 10.1039/D2CC02519A (Communication) Chem. Commun., 2022, 58, 6902-6905

Thermal properties of single-layer MoS2–WS2 alloys enabled by machine-learned interatomic potentials†

Juan M. Marmolejo-Tejada and Martín A. Mosquera *
Department of Chemistry and Biochemistry, Montana State University, Bozeman, MT 59717, USA. E-mail: martinmosquera@montana.edu

Received 4th May 2022 , Accepted 20th May 2022

First published on 20th May 2022

Abstract

Two-dimensional (2D) quantum materials are poised to transform conventional electronics for a wide spectrum of applications that will encompass chemical sciences. For the study of thermal transport in single-layer (1L) or multi-layer transition metal dichalcogenides (TMDs), this work explores the combination of density functional theory (DFT) and algorithmic training for the generation of a moment tensor potential (MTP) that models 1L-MoS2, 1L-WS2 and their alloys, and demonstrates a synergy of theoretical techniques that is anticipated to play an important role in the field. From a high-performance computing perspective, these yield very convenient inter-atomic (or inter-molecular in other contexts) potentials that are useful to predict the response of quantum materials to thermal perturbations, or other driving forces. We show that our trained MTP functions successfully describe vibrational properties of the systems, and their thermal conductivities. The trained potential displays consistent agreement with DFT calculations, as well as the Stillinger–Weber (SW) potential. We also find that the thermal conductivity of the 2D alloys is little affected by sulfur vacancies. This is a behavior that may aid the fine-tuning of material's thermal properties for heat management and energy storage and conversion applications.

Transition metal dichalcogenides (TMDs) are layered two-dimensional (2D) crystalline structures with outstanding electronic and optical properties that have been identified as promising materials for the development of novel nanoelectronic devices, for which accurate estimations of materials' thermal transport properties are essential for thermal management and energy storage/conversion applications.1–9 Structural defects, such as vacancies, usually play a significant role in the electron and thermal transport properties of these devices, which could result in significant lattice distortions in the vicinity of the defect and hinder phonon transport properties.9–11 These systems usually require the preparation of atomistic models with a very large number of atoms (on the order of thousands, or more), for which first-principles methods, such as density functional theory (DFT), become prohibitively expensive in terms of computational cost.

Recently, machine-learned interatomic potentials (MLIPs) have been explored for modeling large-scale molecular systems with high accuracy and significantly lower computational cost than first-principles methods.12,13 Although these potentials usually have very low transferability, they have a flexible functional form and can be systematically improved by increasing the training set and number of basis functions, but at the expense of higher computational cost.12,14–17 Moment tensor potentials (MTPs) describe atomic local environments throughout rotationally covariant tensors (formally, the moment of inertia)12 and have been shown to provide high accuracy and performance, when compared to competing MLIP models in the literature,17 and are expected to help understand the effects of structural and substitutional defects in the performance of TMD-based nanoelectronic devices.13,18

Non-equilibrium molecular dynamics (NEMD) methods have been used for estimating the thermal conductivity, κ, of semiconductor materials due to their direct analogy to experimental procedures.10 For this purpose, the central/scattering region is connected to two semi-infinite, defect-free regions, labeled as the heat source and drain terminals, respectively. Then, the momentum of the hottest atom in the drain is exchanged with the momentum of the coldest atom in the source, although conserving the total energy and momentum of the system.20,21 This produces a heat current from source to drain that results in a temperature gradient along the length of the system.

Alternatively, the non-equilibrium Greens function (NEGF) method has also been employed for calculating thermal conductance, σ, by obtaining the phonon transmission spectrum from the phonon Green's function G(ω)22 and calculating κ = σl/S, where l is the length of the sample and S is the cross-sectional area.23

This work aims at getting a better understanding of the effects of structural and substitutional defects in the performance of nanoelectronic devices and may also provide a tool for phononic engineering by controllably manipulating the thermal and electrical conductivity properties of nanostructured devices,24 where defect engineering may become a necessary step in finely tuning the thermoelectric properties of these devices.8,10,25

Similar to Mortazavi et al.,6,11 we passively train our potential, but include the pristine and alloyed configurations in the training set. This consists of 780 1L configurations from small random atomic displacements and molecular dynamics (MD) trajectories, as detailed in our Supporting Information (ESI†).

To assess the accuracy of our trained MTP, we compute vibrational modes and phonon band diagrams with the supercell method.26–28 For this, we diagonalize the dynamical matrix, which is constructed by small point displacements to the relaxed geometries. We use 21 × 21 × 1 repetitions for 1L-MoS2 and WS2, and 11 × 11 × 1 repetitions for the 1L-MoxW1−xS2 alloys. Then, we calculate the finite differences of the forces with the central finite difference method, using 0.01 Å atomic displacements and 1 × 10−8 Hartree Bohr−2 force tolerance, as implemented in QuantumATK.29,30

Fig. 1 shows the computed phonon band structure diagrams for the 1L-MoxW1−xS2 systems. We find very good agreement between DFT calculations and our MTP. This includes capturing a “negative” energy band at the Γ point in both DFT and MTP calculations; the negative band results from benign computational artifacts and does not translate into any anomalous instability.9,11,31,32 We verify this by fully relaxing the 1L-MoS2 and WS2 unit cells (see Fig. S2, ESI†), where the imaginary frequencies disappear, also noting that the frequency of the optical modes is reduced by ∼20–30 cm−1 with the functional and self-consistent criteria that were used in this study. Interestingly, our trained MTP agrees well with the lattice parameters of the relaxed structures and the resulting phonon band diagrams. To further improve the accuracy of calculations one can add MD trajectories at higher temperatures, relax the unit cells for stress, and increase the accuracy of the reference DFT calculations. However, the random displacement method is sufficient for the modeling performed in this work.


image file: d2cc02519a-f1.tif
Fig. 1 Validation of the trained MTPs: Comparison of different potentials for the determination of phonon dispersion diagrams for 1L- (a) MoS2, (b) WS2, (c) Mo0.75W0.25S2, (d) Mo0.5W0.5S2 and € Mo0.25W0.75S2. These diagrams are calculated with our trained MTP (green-dotted lines) and compared to DFT calculations (solid blue lines). Panel (a) also includes the phonon dispersion of 1L-MoS2 obtained with the SW potential in ref. 19 (orange-dashed lines).

For additional reference, we also compute the phonon dispersion diagram of 1L-MoS2 with the Stillinger–Weber (SW) potential in ref. 19 (see Fig. 1(a)) and compare the optical vibrational modes in Fig. 2(a), finding good agreement with our trained MTP. Furthermore, we verify the evolution of the optical vibrational modes of the compositional-dependent MoxW1−xS2 systems in Fig. 2(b). Consistent with our previous work in ref. 33, we observe a blue-shift of the A2′′ mode (at 477 and 454 cm−1 for pristine 1L-MoS2 and WS2, respectively) and red-shift of the E′′ mode (at 288 and 307 cm−1 for pristine 1L-MoS2 and WS2, respectively) with increased Mo content.


image file: d2cc02519a-f2.tif
Fig. 2 (a) Optical vibrational modes of 1L-MoS2 obtained with our trained MTP and SW potential in ref. 19, and compared to DFT calculations. (b) Evolution of the optical absorption spectra for the 1L-MoxW1−xS2 structures, estimated with our trained MTP (symbols and solid lines) and compared to DFT calculations (dotted lines).

Our inter-atomic potential was then used for estimating the thermal conductivity of 1L-MoS2, WS2 and their alloys, for which we consider sharp and alloy interfaces, as depicted in Fig. 3(a). We setup system configurations of length l = 16, 33, 66 and 132 nm for all cases. For S, we used a fixed width of 25.3 nm for NEMD and 7.58 nm for NEGF calculations (due to larger computational complexity) and assumed a monolayer thickness of 0.6033 nm.34


image file: d2cc02519a-f3.tif
Fig. 3 (a) Top view of the 1L-MoS2/WS2 heterostructure with a sharp or alloy interface. Blue, green and yellow spheres correspond to W, Mo and S atoms, respectively. (b) Phonon thermal conductivity of 1L-MoS2 at T = 300 K as a function of length, using the NEMD and NEGF methods. Our calculations with the trained MTP (solid lines) are compared to calculations with the SW potential (dashed lines) in ref. 19. (c) Phonon thermal conductivity of 1L-MoxW1−xS2 systems at T = 300 K as a function of length, obtained with our trained MTP and the NEMD method. The infinite thermal conductivity, κ∞, in panels (b) and (c) is estimated from the finite length results, assuming a phonon mean free path Λ = 100 nm for all systems. (d) Effect of sulfur vacancies on the thermal conductivity of 1L-MoS2, WS2 and Mo0.5W0.5S2 systems. Alloy interfaces show little variations in the estimated κ with both the NEMD (solid lines) and NEGF (dashed lines) methods.

The conventional approach for estimating the thermal conductivity of infinite systems, κ∞ consists in extrapolating the thermal conductivity results of finite lengths, 1/κ = (1 + Λ/l)/κ∞,1,6,35 where l is the length of the sample and Λ is the effective phonon mean free path, assumed to be 100 nm in all cases.

As summarized in Table S3 (ESI†), previous measurements of κ for 1L-MoS2 at 300 K range between 34.5 ± 4 and 84 ± 17 W K−1m−1,36–38 while theoretical values range between 1.35 and 131 W K−1 m−1, where several methods have been employed, including the self-consistent solution of the Boltzmann transport equation (BTE),28,34,39 and the solution of Slack expression for the acoustic phonon modes,40 in addition to the NEMD41–43 and NEGF23 methods. Our MTP calculations result in κ of 43.2 and 63.4 W K−1 m−1 with the NEMD and NEGF methods, respectively, showing very good agreement with the range of experimental values and our analogous calculations with the SW potential. Our results of finite-length calculations and extrapolation of infinite-length κ are shown in Fig. 3(b).

Fig. 3(c) compares κ for the 1L-MoxW1−xS2 systems using the NEMD method. For 1L-WS2, there is only one reference experimental result for κ to the best of our knowledge, 32 W K−1 m−1,44 which is lower than the experimental references for MoS2. Our MTP calculations result in κ of 33.7 and 51.4 W K−1 m−1 with the NEMD and NEGF methods, respectively, noting that the NEMD result is in better agreement with the experimental value and reference calculations with the Slack expression.40 Then, we verify that 1L-MoS2 has higher κ than 1L-WS2, which is also related to its lower average atomic mass.40

Furthermore, we observe that κ for sharp interfaces is lower than that of the pristine materials, where Mo0.5W0.5S2 has a lower value among the studied concentrations. Alloy interfaces have slightly lower κ and have very low variability among the cases considered in this study, as shown in the inset of Fig. 3(c). This could be useful for engineering the thermal transport properties of other alloy configurations, independent of content concentrations.

Lastly, we add S vacancies to our systems and compute κ for small defect concentrations and a fixed l = 66 nm. As observed in Fig. 3(d), κ is reduced between 25–45% for the pristine materials with a defect density of 0.5%. This is consistent with previous studies where it was observed that vacancies in few-layer MoS2 produced a sharp decrease in κ, that drops even further with increasing defect concentrations due to phonon-defect scattering,25 but noting that very large defect densities result in a transition from a crystalline to an amorphous structure.25 On the other hand, alloy structures are more resilient to S vacancies, with as much as 13% reduction in κ at 0.5% defect density.

We have shown that MTPs can accurately reproduce phononic structure and thermal transport properties of pristine and alloyed 1L-TMDs. Although we find distinctive phonon spectrum features for 1L-MoxW1−xS2 alloys at different concentrations, we observe that thermal transport properties have very low dependence on the alloy content, making it possible to engineer the thermal transport properties of similar materials by mixing different atomic species. Furthermore, we found that alloy configurations are more resilient to S vacancies, providing new opportunities for the fabrication of novel devices without much dependence on costly techniques for obtaining defect-free crystalline materials. Lastly, precision largely depends on the reference DFT method and training set, making it possible to obtain reasonable results without the need for long optimization sequences and higher computational costs.6,11 Future work includes the application of MTPs to multilayer systems and other related quantum materials.

The authors thank Montana State University, Bozeman, for startup support and computational resources within the Tempest Research Cluster, and the Synopsys' QuantumATK development team for helpful discussions. M.A.M. acknowledges support from the MonArk NSF Quantum Foundry supported by the National Science Foundation Q-AMASE-i program under NSF award No. DMR-1906383.

Conflicts of interest

There are no conflicts to declare.

References

  1. P. K. Schelling, S. R. Phillpot and P. Keblinski, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 144306 CrossRef .
  2. J. D. Renteria, D. L. Nika and A. A. Balandin, Appl. Sci., 2014, 4, 525–547 CrossRef .
  3. S. Sun, M. K. Samani, Y. Fu, T. Xu, L. Ye, M. Satwara, K. Jeppson, T. Nilsson, L. Sun and J. Liu, Adv. Mater. Interfaces, 2018, 5, 1800318 CrossRef .
  4. A. Nylander, Y. Fu, M. Huang and J. Liu, IEEE Trans. Compon., Packag., Manuf. Technol., 2019, 9, 427–433 CAS .
  5. Y. Fu, J. Hansson, Y. Liu, S. Chen, A. Zehri, M. K. Samani, N. Wang, Y. Ni, Y. Zhang and Z.-B. Zhang, et al. , 2D Mater., 2019, 7, 012001 CrossRef .
  6. B. Mortazavi, E. V. Podryabinkin, I. S. Novikov, S. Roche, T. Rabczuk, X. Zhuang and A. V. Shapeev, J. Phys.: Mater., 2020, 3, 02LT02 CAS .
  7. G. Bouzerar, S. Thébaud, S. Pecorario and C. Adessi, J. Phys.: Condens. Matter, 2020, 32, 295702 CrossRef CAS PubMed .
  8. A. J. Islam, M. S. Islam, N. Ferdous, J. Park and A. Hashimoto, Phys. Chem. Chem. Phys., 2020, 22, 13592–13602 RSC .
  9. H. Zhan, Y. Nie, Y. Chen, J. M. Bell and Y. Gu, Adv. Funct. Mater., 2020, 30, 1903841 CrossRef CAS .
  10. Y. Lee, S. Lee and G. S. Hwang, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 125202 CrossRef .
  11. B. Mortazavi, I. S. Novikov, E. V. Podryabinkin, S. Roche, T. Rabczuk, A. V. Shapeev and X. Zhuang, Appl. Mater. Today, 2020, 20, 100685 CrossRef .
  12. A. V. Shapeev, Multiscale Model. Simul., 2016, 14, 1153–1173 CrossRef .
  13. V. V. Ladygin, P. Y. Korotaev, A. V. Yanilkin and A. V. Shapeev, Comput. Mater. Sci., 2020, 172, 109333 CrossRef CAS .
  14. J. Behler, J. Chem. Phys., 2016, 145, 170901 CrossRef PubMed .
  15. I. I. Novoselov, A. V. Yanilkin, A. V. Shapeev and E. V. Podryabinkin, Comput. Mater. Sci., 2019, 164, 46–56 CrossRef CAS .
  16. C. Nyshadham, M. Rupp, B. Bekker, A. V. Shapeev, T. Mueller, C. W. Rosenbrock, G. Csányi, D. W. Wingate and G. L.-W. Hart, npj Comput. Mater., 2019, 5, 51 CrossRef .
  17. Y. Zuo, C. Chen, X. Li, Z. Deng, Y. Chen, J. Behler, G. Csányi, A. V. Shapeev, A. P. Thompson, M. A. Wood and S. P. Ong, J. Phys. Chem. A, 2020, 124, 731–745 CrossRef CAS PubMed .
  18. I. S. Novikov, K. Gubaev, E. V. Podryabinkin and A. V. Shapeev, Mach. Learn.: Sci. Technol., 2021, 2, 025002 Search PubMed .
  19. J.-W. Jiang, H. S. Park and T. Rabczuk, J. Appl. Phys., 2013, 114, 064307 CrossRef .
  20. F. Müller-Plathe, J. Chem. Phys., 1997, 106, 6082–6085 CrossRef .
  21. C. Nieto-Draghi and J. B. Avalos, Mol. Phys., 2003, 101, 2303–2307 CrossRef CAS .
  22. T. Markussen, A.-P. Jauho and M. Brandbyge, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 035415 CrossRef .
  23. Y. Cai, J. Lan, G. Zhang and Y.-W. Zhang, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 035438 CrossRef .
  24. G. Xie, Y. Shen, X. Wei, L. Yang, H. Xiao, J. Zhong and G. Zhang, Sci. Rep., 2014, 4, 1–6 Search PubMed .
  25. A. Aiyiti, S. Hu, C. Wang, Q. Xi, Z. Cheng, M. Xia, Y. Ma, J. Wu, J. Guo and Q. Wang, et al. , Nanoscale, 2018, 10, 2727–2734 RSC .
  26. D. Alfè, Comput. Phys. Commun., 2009, 180, 2622–2633 CrossRef .
  27. A. Togo, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 134106 CrossRef .
  28. A. N. Gandi and U. Schwingenschlögl, EPL, 2016, 113, 36002 CrossRef .
  29. QuantumATK version S-2021.06-SP1, Synopsys QuantumATK (https://www.synopsys.com/silicon/quantumatk.html).
  30. S. Smidstrup, T. Markussen, P. Vancraeyveld, J. Wellendorff, J. Schneider, T. Gunst, B. Verstichel, D. Stradi, P. A. Khomyakov and U. G. Vej-Hansen, et al. , J. Phys.: Condens. Matter, 2019, 32(1), 015901 CrossRef PubMed .
  31. V. N. Popov and P. Lambin, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 075427 CrossRef .
  32. H. Zhan, G. Zhang, Y. Zhang, V. Tan, J. M. Bell and Y. Gu, Carbon, 2016, 98, 232–237 CrossRef CAS .
  33. J. M. Marmolejo-Tejada, J. Fix, P. Kung, N. J. Borys and M. A. Mosquera, J. Phys. Chem. C, 2022 DOI:10.1021/acs.jpcc.2c01535 .
  34. W. Li, J. Carrete and N. Mingo, Appl. Phys. Lett., 2013, 103, 253103 CrossRef .
  35. X. Zhang, H. Xie, M. Hu, H. Bao, S. Yue, G. Qin and G. Su, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 054310 CrossRef .
  36. R. Yan, J. R. Simpson, S. Bertolazzi, J. Brivio, M. Watson, X. Wu, A. Kis, T. Luo, A. R. Hight Walker and H. G. Xing, ACS Nano, 2014, 8, 986–993 CrossRef CAS PubMed .
  37. T. Kim, D. Ding, J.-H. Yim, Y.-D. Jho and A. J. Minnich, APL Mater., 2017, 5, 086105 CrossRef .
  38. X. Zhang, D. Sun, Y. Li, G.-H. Lee, X. Cui, D. Chenet, Y. You, T. F. Heinz and J. C. Hone, ACS Appl. Mater. Interfaces, 2015, 7, 25923–25929 CrossRef CAS PubMed .
  39. X. Gu and R. Yang, Appl. Phys. Lett., 2014, 105, 131903 CrossRef .
  40. B. Peng, H. Zhang, H. Shao, Y. Xu, X. Zhang and H. Zhu, RSC Adv., 2016, 6, 5767–5773 RSC .
  41. X. Liu, G. Zhang, Q.-X. Pei and Y.-W. Zhang, Appl. Phys. Lett., 2013, 103, 133113 CrossRef .
  42. V. Varshney, S. S. Patnaik, C. Muratore, A. K. Roy, A. A. Voevodin and B. L. Farmer, Comput. Mater. Sci., 2010, 48, 101–108 CrossRef CAS .
  43. A. Krishnamoorthy, P. Rajak, P. Norouzzadeh, D. J. Singh, R. K. Kalia, A. Nakano and P. Vashishta, AIP Adv., 2019, 9, 035042 CrossRef .
  44. N. Peimyoo, J. Shang, W. Yang, Y. Wang, C. Cong and T. Yu, Nano Res., 2015, 8, 1210–1221 CrossRef CAS .

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cc02519a
This journal is © The Royal Society of Chemistry 2022
Click here to see how this site uses Cookies. View our privacy policy here.