Molecular dynamics investigation on the interfacial thermal resistance between annealed pyrolytic graphite and copper

Abstract

Modern highly integrated microelectronic products often face the challenge of internal heat dissipation, leading to a significant decrease in their operational efficiency. Annealed Pyrolytic Graphite (APG), due to its superior thermal conductivity, has garnered attention from researchers. The interface thermal resistance between APG and supporting materials like copper significantly affects heat transfer during APG's operation. Existing studies rarely delve into the influence of factors such as the shape of APG material interfaces on thermal resistance from a microscopic perspective. In this paper, utilizing transient thermo-reflectance method and non-equilibrium molecular dynamics simulations, the interface thermal resistance of the APG–Cu structure was investigated under different conditions. The impact of parameters such as copper thickness, interface micro-surface morphology, and APG thickness on the calculated interface thermal resistance was examined. Simulation results revealed that copper thickness had a minor effect on the interface thermal resistance. This is because the phonon participation ratio remains unaffected by changes in the thickness of the copper layer. The interfacial thermal resistance beneath microscopically cylindrical copper surfaces was considerably lower than that of rectangular copper surfaces. This is because beneath the cylindrical surface, the enlarged interface contact area facilitates enhanced thermal transport between the interfaces. The computed results of the radial distribution function in the paper also indirectly validate this viewpoint. The magnitude of interfacial thermal resistance for different APG layers was influenced by the coupling effect of intermolecular forces and the layered stacking structure of APG. The interfacial thermal resistance under the condition of three layers of APG reaches its minimum value, which is 2.2 × 10⁻⁹ (K m² W⁻¹). Furthermore, from the phonon perspective, it is found that the interfacial thermal resistance with different numbers of APG layers is closely related to the localization or delocalization state of phonons. As the number of APG layers increased, the interface thermal resistance showed a trend of initial decrease followed by an increase, this is because the average phonon participation ratio increases and then decreases with the number of APG layers. The average phonon participation ratio reaches its maximum value of 0.45778 under the condition of three layers of APG.

---

1. Introduction

In recent years, high-performance electronic products have been evolving towards increased integration.¹ However, this development also leads to a significant accumulation of heat within electronic products, resulting in low operational efficiency.^2–4 Vapor Chambers (VC) is an important method of solving the problem of heat dissipation inside the integrated electronic chip. However, the size and weight of traditional VC are too large, and the internal phase change heat dissipation is prone to induce burn-dry failure and other problems.⁵ Annealed pyrolytic graphite (APG) is an oriented multilayer graphite structure formed by graphene layers through high-temperature chemical vapour deposition process and heat treatment processing above 3000 °C.⁶ APG in plane thermal conductivity can reach up to 1600–2200 W m⁻¹ K⁻¹,⁷ which is four times that of copper (Cu). At the same time, its weight is only a quarter of the weight of Cu, which makes it an ideal material for non-phase-change vapor chambers. However, due to its soft texture, APG needs to be wrapped with metallic materials such as copper to enhance its structural stability. Gurpinar et al.⁸ conducted a study on the heat dissipation of high-power semiconductor modules by using APG instead of copper and found that it can effectively reduce the thermal resistance between the heat source and the cold plate by 50%. Chen et al.⁹ quantitatively analysed the heat dissipation effect of APG/aluminium alloy in the integrated components. Their experimental and simulation results showed that the heat dissipation performance of APG material is much better than that of copper and aluminium alloy plate alone. Su et al.¹⁰ quantitatively analysed the thermal effect of APG in radar signal processing boards and found that, at a total power dissipation of 25 W, the maximum temperature of the aluminium alloy cold plate with APG is reduced by nearly 4 °C and the homogeneity of the temperature is improved from 10 °C to 0.85 °C compared to that of the aluminium alloy cold plate without APG.

APG material enables the transfer of heat through solid-to-solid surface contact, the interfacial thermal resistance between different material interfaces, especially between APG and supporting materials such as copper, has a non-negligible impact on heat transfer.^11,12 Few previous studies have investigated the impact of factors such as the shape of asperities between APG and the substrate material on interfacial thermal resistance from a microscopic point of view. Meanwhile transient thermo-reflectance method, a commonly employed method for assessing the thermal properties of micro- and nano-films,^13,14 is introduced in this study. In the method, two lasers will be used. One induces an instantaneous temperature increase by exposing a small area of the sample surface to a rapidly pulsed laser. The change in optical reflectivity of the sample surface, which is directly related to the surface temperature, is measured using another probe laser beam. From this, the interfacial thermal resistance of the sample can be deduced by analysing the temporal relationship between the laser heating and the temperature response of the sample surface. The advantage of this technique is its high sensitivity to the measurement of the interfacial thermal resistance of small samples or thin films, and its ability to respond well to transient changes in heat. The thermal conductivity between HOPG (Highly Ordered Pyrolytic Graphite), which also has excellent thermal conductivity, and metallic materials has been previously studied by researchers using the transient thermo-reflectance method.¹⁵

In this paper, we intend to use non-equilibrium molecular dynamics (NEMD) simulations to calculate the magnitude of interfacial thermal resistance between APG and the support material Cu based on the transient thermo-reflectance method, i.e., fast pump-probe laser approach.¹⁶ In contrast to the traditional NEMD approach, this pump-probe technique prioritizes examining the dynamic thermal response of the hybrid system, offering significant reductions in computation time. Moreover, the transient thermo-reflectance method has found extensive application in experimental investigations concerning thermal transport in both bulk materials and thin films. The effects of the parameters such as the thickness of the Cu, the microscopic morphology of the interface, and the thickness of the APG on the interfacial thermal resistance calculation will be discussed. Meanwhile, this paper lies in analyzing the influence of copper thickness and the number of APG layers on interfacial thermal resistance from the perspective of phonons. It is also noteworthy that the multi-layer APG conditions in this paper are closer to real-life scenarios. This work aims to provide theoretical data support for the use of APG materials to solve the problems of inefficient chip cooling in highly integrated electronic products, and to provide technical guidance for the design of their material configurations.

2. Mathematical model

2.1 Physical model

In practical applications, APG material is typically encapsulated by the base material Cu. To simplify the model and save computational resources, in this paper, a single-layer APG–Cu interfacial thermal resistance model was firstly developed using Materials Studio,¹⁷ as shown in Fig. 1(a). To investigate the effect of interface microscopic roughness on the interfacial thermal resistance between APG and Cu, two kinds of microscopic surfaces, i.e., cylindrical surface as well as rectangle surface, are constructed for the Cu contact surface in this paper, as shown in Fig. 1(b and c), respectively. Where the parameters d and δ characterise the distance of asperities and depth, respectively. In this paper, d = 2 nm is chosen for the calculation, and the depth δ is chosen as 0.63, 0.83, 1.04, 1.25 and 1.46 nm, respectively. In addition, to be in accordance with the real working conditions, this paper constructs a multilayer APG–Cu interfacial thermal resistance model based on a single-layer APG–Cu interfacial thermal resistance model as shown in Fig. 1(d). We investigate the effect of the number of APG layers on the interfacial thermal resistance of the APG and the Cu.

Fig. 1 APG–Cu interfacial thermal resistance model (a) monolayer APG model (b) cylindrical surface of copper (c) rectangular surface of copper (d) multilayer APG–Cu interfacial thermal resistance model.

The dimensions of the models in this paper: APG: X = 4.2 nm, Y = 18.5 nm, Cu: X = 5.7 nm, Y = 20.1 nm, and Z = 3.8, 4.8, and 5.8 nm, respectively. The lattice length of APG is 0.245999 nm with an angle of 120°, while for Cu, the lattice length is 0.36147 nm with an angle of 90°. The unit cell models mentioned above are all obtained from the Material Studio's model library and have not been artificially modified. In order to ensure the stability of the model and avoid the overlap between the APG layer and the Cu side which may lead to the simulation error report, the distance between the APG and the Cu is set to be 0.346 nm based on the diameter of the C atoms. In the simulation, the X and Y axes of the model are set as periodic boundary conditions, and the Z-axis is set as a free boundary condition. In order to ensure the stability of the model, the cut off radius is set to 1.08 nm in the single-layer APG–Cu interfacial thermal resistance model,¹⁸ and the cut off radius is set to 3.08 nm in the simulation of the multilayer APG–Cu interfacial thermal resistance in order to maximise the reproduction of the real stresses of the multilayer APG–Cu structure.

2.2 Potential function

In this paper, the Embedded-Atom Method (EAM) potential function^19,20 is used to describe the forces between Cu atoms. The EAM potential function expression²¹ is shown in eqn (1):where the EAM potential function considers the effect of parameters such as the force V between particles, the embedding energy F and the local electron density ρ. In addition, N represents the number of atoms in the system and r_ij denotes the distance between two particles i and j. The local electron density expression is shown in eqn (2):where φ represents the correlation term for the local electron density.

In addition, the Lennard-Jones (L-J) potential function¹⁸ is used to describe the forces between C and Cu atoms. The expression for the L–J potential function is shown in eqn (3):

where σ_ij and ε_ij denote the parameters related to distance and energy, respectively, and r_ij is the distance between particles i and j. The parameters of the L–J potential energy between C–Cu atoms are ε = 25.78 meV and σ = 0.30825 nm, respectively.²²

In non-equilibrium molecular dynamics simulations, the Reactive Empirical Bond-order (REBO) potential function, which is adapted from the Tersoff potential function, is usually used to describe the forces between the C atoms inside the APG²³ because the REBO potential function is able to accurately predict the thermophysical properties of graphene as well as the phonon scattering behaviour in the direction perpendicular to the surface of the C atoms (i.e., in the direction of the Z-axis). Its expression is shown in eqn (4):

where E^REBO_ij represents the REBO energy.

2.3 Simulation process

In a transient thermo-reflectance experiment, a laser pulse (pump) is focused on a small spot on the surface of a material. Absorption of the pulse energy by the top material leads to a rapid increase in its temperature, which is then transferred to the substrate by heat conduction, after which the temperature of the material gradually decreases over time. At the same time, the change in temperature of the material leads to a slight change in its optical reflectivity, which can be measured by a second laser pulse (probe). The cooling curve of the material is then obtained and used to determine the interfacial thermal resistance.

In this paper, the molecular dynamics simulations are performed by the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).²⁴ The visualization of the molecular dynamics simulation results is done by Visual Molecular Dynamics (VMD). In the molecular dynamics simulation, the interfacial thermal resistance was calculated as follows: at the beginning of the simulation, the whole system of APG–Cu was relaxed under the canonical ensemble for 1 ns and micro-canonical ensemble for 1 ns, which puts the whole system in an equilibrium state. Time step is 0.5 fs. The temperature-controlling method is Nosé–Hoover method. Subsequently, a heat flow is rapidly introduced into the APG layer, which is equivalent to the first pump laser used to heat the material in the experiment. At the same time, the material temperature variation can be obtained through the simulation, as the probe laser did in real experiment. The temperature and energy of the APG layer as well as the temperature of the Cu side are calculated every 2000 timesteps. The interfacial thermal resistance is calculated by the linear relationship between the energy of the APG layer and the temperature difference at the interface, as shown in eqn (5):

where E_APG represents the energy of the APG layer at different moments, E₀ represents the initial energy of the APG layer at the end of heat flow, A represents the area of the APG layer with the size of 77.7 nm², R represents the interfacial thermal resistance, and the T_APG and T_Cu distributions represent the magnitude of the temperature of the APG layer and the Cu side at different moments.

3. Results and discussion

3.1 Effect of Cu thickness on the interfacial thermal resistance

Matching the magnitude of heat flow in non-metallic materials by transient thermo-reflectance method in ref. 25, a laser irradiated heat flow of 6.24 × 10⁻⁵ W was fed into the single-layer APG–Cu interfacial thermal resistance model for a duration of 50 fs. The energy as well as the temperature changes are shown in Fig. 2. At the initial moment, T_APG = 580 K, T_Cu = 299 K, and E₀ = −3.475 × 10⁻¹⁵ J. As the heat is transferred from the APG layer to the Cu side, the energy of the APG layer starts to decrease rapidly, and correspondingly its temperature starts to decrease gradually. At the same time, the temperature of the Cu side shows a slight increasing. Finally, after 100 ps of thermal relaxation, the temperatures of the APG layer and the Cu side tend to stabilize. To demonstrate the heat import and model temperature distribution more visually in the simulation. The schematic of the interfacial heat flow from APG to Cu is shown in Fig. 3.

Fig. 2 Energy variation and temperature change by the single-layer APG–Cu interfacial thermal resistance model.

Fig. 3 Schematic diagram of computational domain.

Fig. 4 demonstrates the relationship between the energy change of the APG layer and the temperature difference at the interface between APG and Cu, from which we can find that the energy change of the APG layer and the temperature difference at the APG–Cu interface show a good linear relationship. According to eqn (5), the magnitude of the interfacial thermal resistance of the single-layer APG–Cu model can be determined as 1.39 × 10⁻⁸ (K m² W⁻¹). The results of this paper agree with the thermal resistance of graphene–copper interface¹⁸ calculated by Zhang et al. using molecular dynamics simulation.

Fig. 4 The linear relationship between interface temperature difference and energy variation.

Since the scale of molecular dynamics simulation is of micro–nano size, the model size is usually in the range of 1–10 nm.²⁶ Accordingly, in this simulation, the thickness of the Cu side is much smaller than the thickness of the Cu material in the macroscopic experiment. Therefore, it is necessary to explore the influence of the Cu layer thickness on the interfacial thermal resistance. Fig. 5 shows the effect of Cu layer thickness on the interfacial thermal resistance. The data error is analysed after three independent calculations for each condition. From the figure, it can be seen that as the thickness of the Cu layer increases, the interfacial thermal resistance is almost unaffected. In order to deeply investigate the influence of Cu-layer thickness on the interface thermal resistance, this paper introduces the phonon participation ratio (PPR) to analyse the reason why the interface thermal resistance is almost unaffected from the phonon point of view. The PPR²⁷ is weighted by the phonon density of states (PDOS) to obtain the dimensionless number, and its specific expression is shown below:

where P_n(ω) represents the phonon density of states at n frequency and N is the total number of particles. j represents different particles and the formula for the density of phonon states²⁸ is as follows:In the above equation, υ_j(0) is the average velocity vector of the particle at the initial time, υ_j(t) is the velocity of the particle at moment j, and ω represents the vibrational phonon frequency. According to the magnitude of phonon participation, the phonon can be classified as a delocalized mode or a localized mode. The local state phonons are often considered to be not involved in thermal transport. It is well accepted that when R(ω) > 0.3, the phonon is in the delocalized mode, and vice versa, It is in the localized mode. Meanwhile, according to the phonon participation under different Cu-layer thicknesses in Fig. 6, it can be found that the delocalized mode phonons are highly overlapped under different conditions, which is the fundamental reason why the Cu-layer thicknesses have nearly no effect on the interfacial thermal resistance, further confirming the accuracy of the above results. This also means that the Cu-layer thickness Z = 3.8 nm in the initial model is reasonable.

Fig. 5 The influence of copper thickness on interfacial thermal resistance.

Fig. 6 Phonon participation ratio with varying copper thickness.

3.2 Effects of micromorphology on interfacial thermal resistance

In order to investigate the influence of surface morphology on the thermal resistance of APG–Cu interface, cylindrical surface and rectangular surface were constructed in our simulation to simulate the roughness of the micro-interface, respectively. The magnitude of the heat flux is also 6.2487 × 10⁻¹⁶ J. Fig. 7 illustrates the APG–Cu interfacial morphology with different surface configurations at the end of calculation for d = 2 nm and δ = 1.25 nm. From the figure, it can be seen that at the end of the simulation, under the cylindrical surface, the APG layer fits the Cu surface better and has larger the contact area. Under the rectangular surface, the APG layer is always difficult to fit the Cu surface perfectly. The morphology of the APG layer at the final moment does not change significantly compared to the initial moment in Fig. 1. Meanwhile, the surface morphology at different moments is shown in Fig. 8, from which it can be found that under the rectangular surface, the APG layer is always difficult to fit the Cu-side surface. While under the cylindrical surface, the APG layer has already achieved a better fit with the Cu surface within 1 ns, which will increase the contact area between Cu and APG.

Fig. 7 Different surface micromorphology of APG–Cu interface at the end of simulation with a depth of δ = 1.25 nm (a) rectangular surface front view; (b) cylindrical surface front view; (c) cylindrical surface top view; (d) rectangular surface top view.

Fig. 8 Different surface micromorphology of APG–Cu interface at different computational time with a depth of δ = 1.25 nm.

To further demonstrate the increased contact area between APG and Cu beneath the cylindrical surface, this study calculates the Radial Distribution Function (RDF) of different APG–Cu interfacial thermal resistance models with a depth of δ = 1.25 nm. The radial distribution function refers to the probability of finding neighbouring particles around a central particle within a given space. In Fig. 9, we found that the probability of finding the Cu atoms around the C atoms under the cylindrical surface is much larger than that under the rectangular surface, which further confirms that the contact area of APG and Cu is much larger than that under the rectangular surface.

Fig. 9 The radial distribution functions for different surface micromorphology of APG–Cu with a depth δ of 1.25 nm.

Fig. 10 demonstrates the comparison of the magnitude of interfacial thermal resistance under different surface morphology. From the figure, it can be found that at any depth δ, the thermal resistance of the APG–Cu interface under the rectangular surface is much larger than that under the cylindrical surface. Combined with Fig. 7 and 9, we can say that the cylindrical surface shape with better surface fit increases the contact area of the APG–Cu interface and thus reduces the interfacial thermal resistance and improves the heat flow transport efficiency.

Fig. 10 Variation of interfacial thermal resistance as a function of depth δ.

3.3 Effect of APG thickness on interfacial thermal resistance

From a macroscopic point of view, the APG structure cannot be just a two-dimensional structure with a single layer, but a three-dimensional structure with a certain thickness. Thus, it is necessary to investigate the interfacial thermal resistance between multilayer APG and Cu. In this paper, the interfacial thermal resistance between multilayer APG and Cu is simulated and the results are compared with the classical acoustic mismatch model (AMM, acoustic mismatch model) model,²⁹ the H–J model,³⁰ and the diffuse scattering model (DMM, diffuse mismatch model),³¹ as shown in Fig. 11. We can see from the figure that, except the AMM model with poor accuracy, the molecular dynamics simulation in this paper are in good agreement with other classical models, which proves the reliability of the results of the present molecular dynamics simulation.

Fig. 11 Comparison of interfacial thermal resistance between our calculation result and published results.^29–31

In order to further investigate the effect of the number of APG layers on the interfacial thermal resistance, the values of interfacial thermal resistance with 1 to 4 APG layers are simulated, and the results are shown in Fig. 12. From the figure, it can be found that the interfacial thermal resistance under multi-layer APG conditions is smaller than that of single-layer APG conditions. This is attributed to the fact that the molecular force between the interfaces becomes larger as the number of APG layers increases, which enhances the coupling strength between C–Cu and ultimately promotes interfacial heat transport and reduces the interfacial thermal resistance. This is in line with the previous literature.^32,33 It is also worth noting that the calculated interfacial thermal resistance does not decrease all the time as number of APG layers increases. It decreases first and then increase. This is because not only the inter-interfacial molecular forces affect the interfacial thermal resistance, but also the stacked layer structure of APG affects the interfacial thermal transport. Phonons can propagate freely in the single-layer APG structure. While layer–layer interactions in the case of multilayer stacking affect the propagation of phonons and thus hinder the heat transfer.³⁴ Therefore, the intermolecular forces are coupling effect of intermolecular forces and layered stacking structure. The thermal interface resistance in the model exhibits a trend of initial decrease followed by an increase.

Fig. 12 Interfacial thermal resistance as a function of APG layer numbers.

To further elaborate the mechanism of the influence of the number of APG layers on the interfacial thermal resistance, the phonon participation under different numbers of APG layers is calculated and shown in Fig. 13. From the results in the figure, between 5–35 THz, there are fewer phonons in the delocalized mode under single-layer APG, while there are more phonons in the delocalized mode under multi-layer APG participating in thermal transport, and as a result its interfacial thermal resistance is small. Comparing between multilayer APG cases, it is found that there are more delocalized mode phonons around 60 THz in the 3-layer APG compared to the 2-layer and the 4-layer. Therefore, the interfacial thermal resistance is minimum in the case of 3-layer APG.

Fig. 13 The phonon participation under different APG layer numbers.

The average value of phonon participation under different APG layer is shown in Fig. 14, and it can be found that as the number of APG layers increases, the average value of phonon participation increases and then decreases, which means that the contribution of phonons to the heat transfer also gradually increases and then decreases, which corroborates with the results we showed in Fig. 9, and provides a reasonable explanation for the trend of its thermal resistance variation.

Fig. 14 Average phonon participation with varying APG layer numbers.

4. Conclusions

Annealed pyrolytic graphite (APG) has been used by researchers to solve the problem of heat build-up inside modern high-performance electronic products due to its extraordinary thermal conductivity. The interfacial thermal resistance between APG and its contact surfaces has a hard-to-neglect effect on the heat transfer process. In this paper, the magnitude of thermal resistance at the APG–Cu interface is investigated by molecular dynamics simulation. The effect of Cu side thickness on the interfacial thermal resistance is analysed, and the effect of interfacial micromorphology on the thermal resistance of the APG–Cu interface is investigated. The simulation of the thermal resistance of the multilayer APG–Cu interface is carried out. The simulation results show that the magnitude of the interfacial thermal resistance between APG and Cu is almost unchanged at different Cu thicknesses due to the high degree of overlap of the phonons in delocalized mode; the APG–Cu interfacial thermal resistances under the cylindrical surface configurations are much smaller than that under the rectangular surface configurations, which indicates that the fitting degree of the APG layer to the side surface of the Cu determines the interfacial thermal resistance to a great extent. The magnitude of the interfacial thermal resistance is affected by both the intermolecular force and the stacked structure of APG layers. The interfacial thermal resistance shows a decreasing and then increasing with the increase of the number of APG layers. By further analysing the phonon participation in different APG layers, it can be found that more phonons are in the delocalized mode in the multi-layer APG case compared with the single-layer APG case, which promotes the heat transport and leads to a decrease in the interfacial thermal resistance. For the three-layer APG case, the high-frequency phonon participation is larger than the other cases, and there are more phonons in the delocalized mode, so the interfacial thermal resistance is the smallest in this case.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This study was funded by the National Natural Science Foundation of China (No. 52276084).

Notes and references

1. Z. Feng, Z. Jingzhi and H. Xiulan, Int. J. Heat Mass Transfer, 2023, 214, 124453.
2. H. M. Ali, A. Arshad, M. Jabbal and P. G. Verdin, Int. J. Heat Mass Transfer, 2018, 117, 1199–1204.
3. B. Shang, Y. Ma, R. Hu, C. Yuan, J. Hu and X. Luo, Appl. Therm. Eng., 2017, 118, 593–599.
4. M. Janicki and A. Napieralski, Microelectron. J., 2000, 31, 781–785.
5. G. Chen, Y. Tang, S. Zhang, G. Zhong, Y. Sun and J. Li, J. Mech. Eng., 2022, 58, 197–212.
6. A. Moore, A. R. J. P. Ubbelohde and D. Young, Proc. R. Soc. London, Ser. A, 1964, 280, 153–169.
7. C. A. Heusch, H.-G. Moser and A. Kholodenko, Nucl. Instrum. Methods Phys. Res., Sect. A, 2002, 480, 463–469.
8. E. Gurpinar, B. Ozpineci, J. P. Spires and W. Fan, Analysis and Evaluation of Thermally Annealed Pyrolytic Graphite Heat Spreader for Power Modules, 2020 IEEE Applied Power Electronics Conference and Exposition (APEC), 2020, pp. 2741–2747.
9. C. Yixin, S. Lizheng, W. Kejun, L. Jipeng, Zhejiang Huzhou, Research on Integrated Design and Technology of Annealed Pyrolytic Graphite Structure, The 16th China CAE Annual Conference, 2020, vol. 4, pp. 237–240.
10. Y. X. Chen, L. Z. Su and K. J. Wang, et al., Thermal Conductivity Analysis of APG and Its Application in Radar, Fire Control Radar Technology, 2021, 50, pp. 90–94.
11. W. Shengtao, Z. Kejie and L. Heng, Commun. Technol., 2013, 46, 101–104.
12. D. Guozhong, D. Zechun, L. Guoyou and L. Jinhui, Power Electron. Technol., 2018, 52, 34–37.
13. Y. Xu, H. Wang, Y. Tanaka, M. Shimono and M. Yamazaki, Mater. Trans., 2007, 48, 148–150.
14. R. Kato and I. Hatta, Int. J. Thermophys., 2008, 29, 2062–2071.
15. A. J. Schmidt, K. C. Collins, A. J. Minnich and G. Chen, J. Appl. Phys., 2010, 107, 104907.
16. R. Stoner and H. Maris, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 16373.
17. L. Jiaxuan, Metall. Mater., 2021, 41, 62–63.
18. Y. Hong, L. Li, X. C. Zeng and J. Zhang, Nanoscale, 2015, 7, 6286–6294.
19. M. I. Mendelev, S. Han, D. J. Srolovitz, G. J. Ackland, D. Y. Sun and M. Asta, Philos. Mag., 2003, 83, 3977–3994.
20. G. Bonny, R. C. Pasianot, N. Castin and L. Malerba, Philos. Mag., 2009, 89, 3531–3546.
21. M. S. Daw and M. I. Baskes, Phys. Rev. Lett., 1983, 50, 1285–1288.
22. L. Lindsay, D. Broido and N. Mingo, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 1321–1330.
23. J. Tersoff, Phys. Rev. Lett., 1988, 61, 2879.
24. L. Zhao, Research on the Thermophysical Properties Enhancement of Molten Salt Based Nanofluid, North China Electric Power University, Beijing, p. 2021.
25. S. Varghese, J. D. Mehew, A. Block, D. S. Reig, P. Woźniak, R. Farris, Z. Zanolli, P. Ordejón, M. J. Verstraete and N. F. Van Hulst, Rev. Sci. Instrum., 2023, 94, 012005.
26. S. Gowthaman, Funct. Compos. Struct., 2023, 5, 012005.
27. Y. Guo, M. Bescond, Z. Zhang, S. Xiong, K. Hirakawa, M. Nomura and S. Volz, APL Mater., 2021, 9, 091104.
28. Y. Yang, J. Ma, J. Yang and Y. Zhang, ACS Appl. Mater. Interfaces, 2022, 14, 45742–45751.
29. W. Little, Can. J. Phys., 1959, 37, 334–349.
30. H. Cao, D.-B. Xiong, Z. Tan, G. Fan, Z. Li, Q. Guo, Y. Su, C. Guo and D. Zhang, J. Alloys Compd., 2019, 771, 228–237.
31. E. T. Swartz and R. O. Pohl, Rev. Mod. Phys., 1989, 61, 605.
32. Z. Ding, Q.-X. Pei, J.-W. Jiang, W. Huang and Y.-W. Zhang, Carbon, 2016, 96, 888–896.
33. A. I. Khan, R. Paul and S. Subrina, RSC Adv., 2017, 7, 44780–44787.
34. L. Dongjing, Z. Fu, C. Shuaiyang and H. Zhiliang, Acta. Phys. Sin., 2023, 72, 157901.

---