Quasiballistic thermal transport in submicron-scale graphene nanoribbons at room-temperature

Abstract

Phonon transport in two-dimensional materials has been the subject of intensive studies both theoretically and experimentally. Recently observed unique phenomena such as Poiseuille flow at low temperature in graphene nanoribbons (GNRs) initiated strong interest in similar effects at higher temperatures. Here, we carry out massive molecular dynamics simulations to examine thermal transport in GNRs at room temperature (RT) and demonstrate that non-diffusive behaviors including Poiseuille-like local thermal conductivity and second sound are obtained, indicating quasiballistic thermal transport. For narrow GNRs, a Poiseuille-like thermal conductivity profile develops across the nanoribbon width, and wider GNRs exhibit a mixed nature of phonon transport in that diffusive transport is dominant in the middle region whereas non-uniform behavior is observed near lateral GNR boundaries. In addition, transient heating simulations reveal that the driftless second sound can propagate through GNRs regardless of the GNR width. By decomposing the atomic motion into out-of-plane and in-plane modes, it is further shown that the observed quasiballistic thermal transport is primarily contributed by the out-of-plane motion of C atoms in GNRs.

---

Introduction

Advances in nanoscience and nanotechnology have enabled thermal conductivity (κ) to be studied from a geometry-dependent viewpoint rather than taking κ to be an intrinsic property as is done for the bulk. For instance, graphitic materials such as graphene nanoribbons (GNRs) and carbon nanotubes (CNTs) have displayed size-dependent thermal conductivity which becomes logarithmically divergent with the sample length.^1–7 This abnormality in thermal conductivity of graphene nanostructures has been attributed to the unique phonon transport originating from out-of-plane acoustic (ZA) phonon branches which possess particularly long mean free paths (MFPs).^7–9 Moreover, previous studies on graphene reported that the contribution of ZA phonon branches to thermal conductivity becomes more than 70% of total thermal conductivity even after including four-phonon scattering processes which are known to further reduce thermal conductivity in graphene.^10,11 Such distinctive characteristics in MFPs have indeed made ZA phonons a primary contributor to high thermal conductivity in graphene systems.^12–14

With the important role played by ZA phonons, thermal transport in graphene nanostructures is normally studied from two distinct perspectives. In one regime occurring at low temperatures, phonon MFPs are considerably larger than sample size so that phonons pass through the sample without suffering from significant scattering except at the boundaries, which results in ballistic phonon transport. The other is the high-temperature diffusive regime where all phonon scattering events are included because phonon MFPs are much smaller than sample size. Such a dichotomic approach based on contrasting length scales has been particularly successful such that it has provided valuable understanding and insights into phonon transport in nanostructured systems. Meanwhile, recent years have witnessed high interest in phonon transport occurring between ballistic and diffusive regimes. In these intermediate regimes, which have received considerable attention concerning thermal transport in two-dimensional (2D) materials both theoretically^15–22 and experimentally,^23–26 thermal transport exhibits unique behaviors that are different from ballistic or diffusive regimes. One is hydrodynamic phonon transport occurring at low temperatures,^26–28 in which normal scattering between phonons is highly dominant and phonon momentum is conserved during the scattering processes until it is eventually dissipated by other resistive scattering mechanisms. In the hydrodynamic regime, phonon transport displays two distinctive phenomena in thermal propagation, Poiseuille flow and second sound.^27,29–31 The Poiseuille flow is a non-uniform heat flux in a steady-state and mainly caused by diffuse scattering at the boundaries, and second sound is a non-Fourier heat conduction where the heat propagates in a wavelike form in transient heating. Indeed, Poiseuille flow and second sound have been observed in diverse, submillimeter-sized samples including graphene nanostructures at low temperatures.^28,32–36

It is noted that Poiseuille-like flow and second sound are also achievable for another type of intermediate scale phonon transport, the quasiballistic regime. In cases of quasiballistic transport, some of the phonons propagate diffusively through multiple scattering events while others are transported ballistically.^37–40 Quasiballistic transport, which occurs at similar or higher temperatures than the hydrodynamic one, concerns the transition between ballistic and diffusive thermal transport, and such crossover is shown to take place when sample sizes are allowed to vary.^40,41 For instance, the total thermal conductivities of Ni lines, Si and Ge were seen to exhibit different contributions from ballistic and diffusive transport. While these studies shed light on thermal transport within the intermediate regime, a detailed understanding of ballistic and diffusive contributions such as the local distribution of the two effects inside samples and roles of distinct vibrational modes in phonon transport is lacking although such information will be highly valuable both fundamentally and technologically.

To address this issue, in the present work we carry out large-scale molecular dynamics (MD) simulations on graphene nanoribbons (GNRs) to investigate their thermal transport behavior at room temperature (RT) by systematically varying the width and length. It is demonstrated from the simulations that phonon transport in GNRs reveals quasiballistic behavior at RT, which is carefully analyzed by examining the local thermal conductivity profile and heat pulse propagation. Moreover, in-plane and out-of-plane vibrational modes are shown to make different contributions to thermal transport in GNRs. These findings surely provide a useful insight into phonon transport in graphene nanoribbons, which may lead to advancements in thermal management applications.

Computational details

All MD simulations are performed with the GPUMD (Graphics Processing Units Molecular Dynamics, version 2.5.1) open source program,⁴² and the Tersoff potential optimized for graphitic C–C covalent systems is used to describe interatomic interaction.⁴³ Periodic boundary conditions are applied to the heat flow (x) direction, whereas both width (y) and thickness (z) directions are treated as free. To examine the size effect on thermal transport in GNRs, we separately vary their width (W) and length (L) ranging from 10 to 1000 nm, which leads to as many as 50 million C atoms within a simulation cell.

Steady-state heating conditions

Fig. 1(a) shows the geometry of GNRs employed in the present work. The entire length of GNRs is divided into 4 subdivisions which are denoted as a fixed, heat source, heat flow and heat sink region, respectively, in the figure. A fixed region, which consists of stationary 5 zigzag C chains across the ribbon width, is introduced to satisfy periodic boundary conditions, and the heat flow region of length L is further divided into 59 subgroups along the length for detailed examination of heat flux. To investigate thermal transport, we carry out nonequilibrium molecular dynamics (NEMD) simulations at 300 K with a time step of 1 fs and employ the Langevin thermostat to generate the temperature difference (ΔT) between the heat source and heat sink. The resulting temperature difference will then bring about the heat flux q which is given asHere, Q(= |dE/dt|) is the energy exchange rate between the groups in GNRs and S is the cross-sectional area of the nanoribbon perpendicular to the heat flow direction, respectively. Once the steady-state is reached with the intended ΔT, the thermal conductivity of a GNR with length L is computed from q as

κ = (q/ΔT)L (2)

In this study, ΔT is set to be 60 K and total production time is as long as 50 ns depending on the system size.

Fig. 1 (a) Structure of graphene nanoribbons for MD simulations. L denotes the length of the heat flow region between the heat source and heat sink, and W is the ribbon width, respectively. (b) Thermal conductivity of GNRs as a function of L for different W values.

Transient heating conditions

In order to examine the propagation of heat pulses, we consider the GNR with (W, L) = (250 nm, 1000 nm) as a representative system and keep track of the temperature profile as a function of time after an initial heat pulse. Specifically, by taking the reference temperature (T_ref) and time step to be 300 K and 0.5 fs, respectively, the heat pulse, which is modeled as a 90% higher temperature than T_ref, is imposed on the heat source region for 0.125 ps considering phonon relaxation time.⁴⁴ In transient heating cases, the length of the heat source region is reduced to that of a single armchair unit of graphene to properly take the localized nature of the heat pulse into account. For systematic analysis of the propagating behavior of the heat pulse, the heat flow region of the GNR is divided into much finer subgroups (2456 subgroups) than steady-state heating cases, and the following temperature profile, which is computed from the equipartition theorem for each group along the length, is monitored along the GNR length within NVE simulations.

Results and discussion

Lattice thermal conductivity under steady-state heating

We first calculate the lattice thermal conductivity of the GNRs. To assess the effect of the nanoribbon size on thermal conductivity, we consider the GNR lengths ranging from 10 to 1000 nm at each of the W = 100, 250, 500 and 1000 nm cases. As is seen from Fig. 1(b), κ shows a logarithmic increase with the length of GNRs as κ ∼ a(log L)^b, where b increases from 1.96 when W = 100 nm and becomes as high as 2.31 for W = 1000 nm. It is noted that the logarithmic and power-law dependence of κ on L are often observed in 2D and 1D systems, respectively.^45–48 Due to the finite sizes in both length and width dimensions, the thermal conductivity of the present systems can be understood by phenomenologically combining logarithmic and power-law behaviors.

To achieve an additional insight into the width dependence of κ, we examine the position-dependent thermal conductivity in the width direction (y direction in Fig. 1(a)) for different W and L cases. Fig. 2 shows the position-dependent thermal conductivity across the width (referred to as “local thermal conductivity” hereafter),²¹κ(ȳ) with ȳ = y/W, which is averaged over the ribbon length for each of (W, L) values. It is first noted from Fig. 2(a)–(c) that regardless of the width, there exist flat regions in κ(ȳ) for the shortest graphene nanoribbons and the flat region becomes more pronounced for wider GNRs. For instance, while the local thermal conductivity of the GNR with (W, L) = (250 nm, 100 nm) tends to show a slightly concave curvature (Fig. 2(a)), κ(ȳ) of the 1000 nm-wide GNR with the same length reveals an almost completely flat profile in most regions across the width (Fig. 2(c)). As the ribbon length grows, however, the flat region in κ(ȳ) becomes narrower and eventually disappears except in the widest cases. When W = 250 nm, the flat behavior nearly disappears even for the 250 nm-long GNR, whereas it remains observable for as long as 500 nm in W = 500 nm cases and persists for all lengths in 1000 nm-wide nanoribbons, as is seen in Fig. 2(b) and (c), respectively. It is noted that this behavior is similar to that of hydrodynamic phonon transport,²⁷ but the present result is differentiated from hydrodynamic transport since thermal transport is additionally contributed by diffusive characters as discussed below. Also, collective phonon excitation (CPE) is known to play a significant role in transport.⁴ It should be noted, however, that CPE is characterized by a hundred-micrometer-long mean free path, so in 1 μm-sized GNRs its contribution will not be significant.

Fig. 2 Length-averaged position-dependent thermal conductivity κ (in units of W mK⁻¹) along width: (a) W = 250 nm, (b) 500 nm and (c) 1000 nm, respectively. The black dashed lines represent the fitting function in eqn (3).

Notably, the computed local thermal conductivity is perfectly fitted with the following function from the phonon hydrodynamic equation for nanoscale heat transport, which is developed to understand quasiballistic behavior²² where λ, Kn, and s are the fitting parameters, respectively.

Here, Kn is the Knudsen number which is defined in terms of phonon mean free path (MFP, Λ) and GNR width as Kn = Λ/W. It is particularly interesting from Fig. 2(a) and (b) that in cases of W = 250 nm and 500 nm, κ(ȳ) shows Poiseuille-flow-like profiles for the longest GNRs, which is seen from the fact that the local thermal conductivity closely follows f(ȳ) = (β − ȳ)^η(β + ȳ)^η with η being as high as 0.62. Moreover, there exist damping effects around ȳ = 0 and 1 in thermal conductivity for all widths and lengths, which arise from the interaction between phonons and GNR boundaries. Interestingly, the damping effects become stronger as the GNR grows longer, which is evidenced by the observation that κ_cb = [small kappa, Greek, macron](0.5)/[small kappa, Greek, macron](0), and the thermal conductivity ratio between ȳ = 0 and 0.5 increases from 1.20 at L = 100 nm to 1.30 at L = 1000 nm. Clearly, all these characteristics in the local thermal conductivity point to the quasiballistic nature of phonon transport in graphene nanoribbons, which is also evidenced by Kn values. From eqn (3), it is found that Kn ranges from 0.1 to 0.6 for most GNRs, which indicates that the corresponding phonon transport lies between ballistic and diffusive regimes.^29,49 Interestingly, a deeper understanding of the quasiballistic nature of the phonon transport in GNRs can be attained by investigating the local thermal conductivity at different locations along the GNR length.

To this end, κ(ȳ) is calculated for each group within the heat flow region of 1000 nm-long GNRs and plotted for every third group as shown in Fig. 3, where open circles represent the Knudsen layer thickness (δ_K). Here, δ_K is set to be Λ,^22,50 which is extracted from the Knudsen number by applying the same fitting function (eqn (3)) to κ(ȳ) in each group. As is seen from Fig. 3(a), the Knudsen layers of the 250 nm-wide GNR rapidly develop from each boundary and eventually overlap about 150 nm away from the heat source, and κ(ȳ) maintains its profile thereafter. This phenomenon is very similar to a fluid flow in a pipe where the velocity of an entering flow shows a fully developed form after a finite distance from the pipe entrance. Unlike ordinary fluids, the Knudsen layer begins to split again around 100 nm from the heat sink and each branch narrows as the heat sink is further approached. However, the wider GNRs reveal different behaviors in the thermal conductivity from the 250 nm-wide cases. For the GNR of 500 nm width, it is seen from Fig. 3(b) that while Knudsen layers again gradually increase as heat is transported from the heat source, they do not completely merge with each other, which is seen by observing that δ_K holds a maximum value halfway along the length and decreases again thereafter. κ(ȳ) of the 1000 nm-wide GNR also shows the same behavior as is presented in Fig. 3(c).

Fig. 3 Local thermal conductivity at each location along the length for (a) 250 nm-, (b) 500 nm- and (c) 1000 nm-wide GNRs, respectively. The length of GNRs is fixed at 1000 nm. Open circles denote the Knudsen layer thickness fitted with black dotted curves. Top panels are magnified figures of thermal conductivity near the heat source and heat sink with the results for more groups included (gray curves), and the bottom panel represents the heat flow region in Fig. 1(a).

It is noted that the Knudsen layer represents a region within which phonon-boundary scattering plays an important role and thus produces a viscous damping effect. Consequently, thermal conductivity is reduced with a spatial gradient and Poiseuille-flow-like behavior thus appears in phonon transport which is similar to gas flow in pipes.^21,29 Rapid merge of the Knudsen layers in the 250 nm-wide GNR in Fig. 3(a) indicates that the phonon transport in this case is predominantly determined by the boundary-induced damping effect. On the other hand, 500 nm- and 1000 nm-wide GNRs exhibit a region which is not influenced by the damping effect as is presented in Fig. 3(b) and (c). The average values of δ_K are 164.0 and 173.0 nm for 500 nm- and 1000 nm-wide GNRs with maxima of 206.0 and 221.4 nm, respectively, resulting in 2δ_K < W for these wider GNRs. Thus, κ(ȳ) always shows a flat profile of a finite size (l_F) across the width in these nanoribbons, which points to diffusive phonon transport governed by Fourier's law.²² This observation demonstrates the transitional behavior in phonon transport of GNRs as evidenced by the coexistence of quasiballistic and diffusive characteristics. Moreover, it is the relative ratio between δ_K and l_F that determines the overall strength of quasiballistic phonon transport in GNRs, for which the structural parameters (W and L) play an important role.

What is interesting is that a further insight into the thermal transport in GNRs can be achieved by decomposing κ(ȳ) into in-plane (κ_i(ȳ)) and out-of-plane (κ_o(ȳ)) components.⁵¹Fig. 4(a)–(c) show the decomposed local thermal conductivity of the longest GNRs for different widths. As is clearly seen from the figures, κ_i(ȳ) is almost flat across the width for all cases whereas κ_o(ȳ) retains nearly the same profile as the entire thermal conductivity. Thus, the in-plane thermal transport becomes mostly diffusive and the quasiballistic nature of the phonon transport in GNRs is governed predominantly by the out-of-plane component of thermal conductivity.

Fig. 4 Local thermal conductivity κ (in units of W mK⁻¹) decomposed into in-plane (κ_i, blue) and out-of-plane (κ_o, red) contributions for different widths of GNRs: (a) W = 250 nm, (b) 500 nm and (c) 1000 nm, respectively. The length of the nanoribbons is fixed at 1000 nm. The black dashed lines represent the fitting function in eqn (3).

Heat waves under transient heating

To further understand the characteristics of phonon transport, we examine the propagation of the heat pulse which is simulated through transient heating. Fig. 5(a) presents the temperature profiles of the GNR with (W, L) = (250 nm, 1000 nm) at T_ref = 300 K along the length until 7 ps after the heat pulse is imposed. As is seen from the figure, the initial pulse at x = 0 begins to propagate towards the heat sink, and more interestingly, there exists a distinct peak in the temperature profile, which first appears at x ∼ 10 nm and then travels at a constant speed with a reduced height. However, the propagation of heat pulses at 900 K shows an apparent difference from that of 300 K cases as is seen from Fig. 5(b). While a peak is initially observed in the temperature profile around the same position as in 300 K, the height is much smaller: it is only 0.8% higher than the surrounding temperature at T_ref = 900 K, which is compared to 2.5% for 300 K. Moreover, the peak rapidly disappears through dissipation into the background as early as 5 ps and the subsequent thermal propagation is governed by Fourier's law for heat conduction, which can be easily seen from the linear decrease in the temperature profile (Fig. 5(b)).

Fig. 5 The propagating temperature profiles along the length in GNRs with (W, L) = (250 nm, 1000 nm) for (a) T_ref = 300 K and (b) 900 K at different times, respectively. The decomposed temperature profiles at different propagating times: (c) out-of-plane and (d) in-plane temperature for T_ref = 300 K, and (e) out-of-plane and (f) in-plane temperature for T_ref = 900 K, respectively.

It is noted that a better understanding of heat propagation is achieved by examining mode-dependent temperatures as is done for steady-state heating cases. To this end, the temperature profile is further decomposed into in-plane (T_i) and out-of-plane (T_o) contributions, which are defined from the equipartition theorem as and , respectively. Here, N and k_B are the total number of atoms in a given group and Boltzmann constant and M_j and v_j,α are the mass and velocity in the α-direction (α = x, y, z) of the j-th atom, respectively. This decomposition is relevant since ZA phonons make a significant contribution to thermal transport.¹¹Fig. 5(c) and (d) show T_o and T_i at 300 K, respectively, and as is shown in the figures, the peak (solid circles in Fig. 5(c)) is mostly contributed by the out-of-plane motion of the atoms in GNRs with a minor addition from the in-plane velocity. What is interesting is the existence of a shoulder structure (SS) in T_i (dotted circles in Fig. 5(d)), which is moving slightly ahead of the peak and persists until 3 ps. Noticeably, both the peak and shoulder structure in T_i do not exist after 3 ps since the thermal transport becomes diffusive. It is found that the propagation speeds of the peak and SS are 8.5 km s⁻¹ and 11.9 km s⁻¹, respectively, which indicates that the former and latter are generated by the second sound and ballistic transport, respectively. Interestingly, a recent study based on the modal energy analysis showed that second sound and ballistic pulse are directly characterized by out-of-plane mode (ZA branch) and in-plane mode (dominantly from the TA branch), respectively.⁵² However, insufficient normal scattering in the present cases causes the obtained second sound to be driftless as in previous studies,^31,52–54 in which the heat pulses in submicron-scale GNRs are mainly carried by the driftless second sound from ZA phonons while the ballistic contributions from in-plane phonon modes clearly suffer from dissipation. At a higher temperature of T_ref = 900 K, however, while the small peak in T_o is still present after the initial temperature pulse, it rapidly dissipates into background after existing only for as long as 3 ps (Fig. 5(e)). Also, as is seen from the T_i profiles in Fig. 5(f), the lifetimes of the peak and SS are around 1 ps, which is much shorter than those for T_ref = 300 K cases. These findings clearly imply that the driftless second sound and ballisticity are hardly existent at 900 K and phonon transport is primarily diffusive.

Fig. 6 shows the decomposed temperature profiles for (W, L) = (500 nm, 500 nm) (Fig. 6(a) and (b)) and (1000 nm, 100 nm) (Fig. 6(c) and (d)), respectively. It is clearly seen from the figures that both the T_i and T_o profiles of these wider GNRs are nearly the same as those for (W, L) = (250 nm, 1000 nm). In particular, the position and height of the peak in T_o are well preserved and the peculiar shoulder structure in T_i is maintained, which implies that the nature of phonon transport in submicrometer-sized GRNs is well preserved at RT regardless of structural variations.

Fig. 6 The decomposed temperature profiles in GNRs along the length at different propagating times: (a) out-of-plane and (b) in-plane temperature for (W, L)=(500 nm, 500 nm), and (c) out-of-plane and (d) in-plane temperature for (W, L)=(1000 nm, 100 nm), respectively. Simulations are conducted at T_ref = 300 K.

Conclusion

We have explored the phonon thermal transport of submicron-scale graphene nanoribbons at room temperature with a range of widths and lengths by employing massive molecular dynamics simulations. The simulation results clearly demonstrate that phonon transport in these GNRs becomes mostly quasiballistic, as is evidenced by the presence of non-diffusive local thermal conductivity and driftless second sound. In particular, non-uniform thermal transport is observable within the Knudsen layer from the lateral boundary of GNRs, which makes the Poiseuille-like profile in local thermal conductivity visible only for narrow ribbons. More importantly, our MD simulations reveal that the out-of-plane motion of C atoms plays a crucial role in generating quasiballistic behavior in thermal transport whereas diffusive transport is contributed by in-plane atomic vibration. These findings will shed new light on high-temperature phonon transport and advance the understanding of different contributions of distinctive vibrating modes to heat conduction in nanostructured materials.

Conflicts of interest

The authors declare no competing interests.

Acknowledgements

This work was supported by the GIST-MIT Collaboration grant funded by GIST in 2024. JHL also acknowledges support from the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korean Government (20212050100010).

References

1. X. Xu, L. F. C. Pereira, Y. Wang, J. Wu, K. Zhang, X. Zhao, S. Bae, C. Tinh Bui, R. Xie, J. T. L. Thong, B. H. Hong, K. P. Loh, D. Donadio, B. Li and B. Özyilmaz, Nat. Commun., 2014, 5, 3689.
2. A. K. Majee and Z. Aksamija, Phys. Rev. B, 2016, 93, 235423.
3. G. Barbalinardo, Z. Chen, H. Dong, Z. Fan and D. Donadio, Phys. Rev. Lett., 2021, 127, 025902.
4. G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari and F. Mauri, Nano Lett., 2014, 14, 6109–6114.
5. Z. Zhang, Y. Ouyang, Y. Guo, T. Nakayama, M. Nomura, S. Volz and J. Chen, Phys. Rev. B, 2020, 102, 195302.
6. Z. Zhang, Y. Ouyang, Y. Cheng, J. Chen, N. Li and G. Zhang, Phys. Rep., 2020, 860, 1–26.
7. D. L. Nika, A. S. Askerov and A. A. Balandin, Nano Lett., 2012, 12, 3238–3244.
8. S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau and A. A. Balandin, Nat. Mater., 2010, 9, 555–558.
9. A. A. Balandin, Nat. Mater., 2011, 10, 569–581.
10. T. Feng and X. Ruan, Phys. Rev. B, 2018, 97, 45202.
11. Z. Han and X. Ruan, Phys. Rev. B, 2023, 108, L121412.
12. L. Lindsay, D. A. Broido and N. Mingo, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 115427.
13. J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, N. Mingo, R. S. Ruoff and L. Shi, Science, 2010, 328, 213–216.
14. D. L. Nika and A. A. Balandin, J. Phys.: Condens. Matter, 2012, 24, 233203.
15. S. Lee, D. Broido, K. Esfarjani and G. Chen, Nat. Commun., 2015, 6, 6290.
16. S. Lee and L. Lindsay, Phys. Rev. B, 2017, 95, 184304.
17. C. Zhang and Z. Guo, Int. J. Heat Mass Transfer, 2021, 181, 121847.
18. X.-P. Luo, Y.-Y. Guo, M.-R. Wang and H.-L. Yi, Phys. Rev. B, 2019, 100, 155401.
19. Y. Guo, Z. Zhang, M. Bescond, S. Xiong, M. Wang, M. Nomura and S. Volz, Phys. Rev. B, 2021, 104, 075450.
20. Y. Guo and M. Wang, Int. J. Therm. Sci., 2022, 171, 107178.
21. X. Li and S. Lee, Phys. Rev. B, 2018, 97, 094309.
22. Y. Guo and M. Wang, Phys. Rev. B, 2018, 97, 035421.
23. S. Huberman, R. A. Duncan, K. Chen, B. Song, V. Chiloyan, Z. Ding, A. A. Maznev, G. Chen and K. A. Nelson, Science, 2019, 364, 375–379.
24. Z. Ding, K. Chen, B. Song, J. Shin, A. A. Maznev, K. A. Nelson and G. Chen, Nat. Commun., 2022, 13, 285.
25. Y. Machida, N. Matsumoto, T. Isono and K. Behnia, Science, 2020, 367, 309–312.
26. Y. Machida, A. Subedi, K. Akiba, A. Miyake, M. Tokunaga, Y. Akahama, K. Izawa and K. Behnia, Sci. Adv., 2023, 4, eaat3374.
27. A. Cepellotti, G. Fugallo, L. Paulatto, M. Lazzeri, F. Mauri and N. Marzari, Nat. Commun., 2015, 6, 6400.
28. X. Huang, Y. Guo, Y. Wu, S. Masubuchi, K. Watanabe, T. Taniguchi, Z. Zhang, S. Volz, T. Machida and M. Nomura, Nat. Commun., 2023, 14, 2044.
29. Y. Guo and M. Wang, Phys. Rep., 2015, 595, 1–44.
30. C. Yu, Y. Ouyang and J. Chen, J. Appl. Phys., 2021, 130, 010902.
31. B.-D. Nie and B.-Y. Cao, Nanoscale Microscale Thermophys. Eng., 2020, 24, 94–122.
32. R. A. Guyer and J. A. Krumhansl, Phys. Rev., 1966, 148, 778–788.
33. X. Li and S. Lee, Phys. Rev. B, 2019, 99, 085202.
34. K. Ghosh, A. Kusiak and J.-L. Battaglia, Phys. Rev. Mater., 2021, 5, 073605.
35. M. Markov, J. Sjakste, G. Barbarino, G. Fugallo, L. Paulatto, M. Lazzeri, F. Mauri and N. Vast, Phys. Rev. Lett., 2018, 120, 075901.
36. A. Beardo, M. López-Suárez, L. A. Pérez, L. Sendra, M. I. Alonso, C. Melis, J. Bafaluy, J. Camacho, L. Colombo, R. Rurali, F. X. Alvarez and J. S. Reparaz, Sci. Adv., 2023, 7, eabg4677.
37. M.-H. Bae, Z. Li, Z. Aksamija, P. N. Martin, F. Xiong, Z.-Y. Ong, I. Knezevic and E. Pop, Nat. Commun., 2013, 4, 1734.
38. A. Sood, F. Xiong, S. Chen, R. Cheaito, F. Lian, M. Asheghi, Y. Cui, D. Donadio, K. E. Goodson and E. Pop, Nano Lett., 2019, 19, 2434–2442.
39. R. Anufriev, S. Gluchko, S. Volz and M. Nomura, ACS Nano, 2018, 12, 11928–11935.
40. P. Torres, M. Royo, M. López-Suárez, J. Shiomi and R. Rurali, Phys. Rev. B, 2020, 102, 144305.
41. M. E. Siemens, Q. Li, R. Yang, K. A. Nelson, E. H. Anderson, M. M. Murnane and H. C. Kapteyn, Nat. Mater., 2010, 9, 26–30.
42. Z. Fan, W. Chen, V. Vierimaa and A. Harju, Comput. Phys. Commun., 2017, 218, 10–16.
43. L. Lindsay and D. A. Broido, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 205441.
44. B. Qiu and X. Ruan, Appl. Phys. Lett., 2012, 100, 193101.
45. L. Wang, B. Hu and B. Li, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 40101.
46. L. Wang and T. Wang, Europhys. Lett., 2011, 93, 54002.
47. Y. Deng and S. W. Cranford, Comput. Mater. Sci., 2017, 129, 226–230.
48. A. G. Remanidevi and D. Barik, Phys. Rev. E, 2019, 99, 022103.
49. A. Sellitto, I. Carlomagno and D. Jou, Proc. R. Soc. A, 2015, 471, 20150376.
50. G. H. Tang, Y. Zhao, G. X. Zhai and C. Bi, J. Appl. Phys., 2011, 110, 046102.
51. Z. Fan, L. F. C. Pereira, P. Hirvonen, M. M. Ervasti, K. R. Elder, D. Donadio, T. Ala-Nissila and A. Harju, Phys. Rev. B, 2017, 95, 144309.
52. C. Yu, S. Shan, S. Lu, Z. Zhang and J. Chen, Phys. Rev. B, 2023, 107, 165424.
53. W.-J. Yao and B.-Y. Cao, Chin. Sci. Bull., 2014, 59, 3495–3503.
54. H. Rezgui, ACS Omega, 2023, 8, 23964–23974.

---