AAPPS bulletin

Research and Review

Phonon thermal transport in bilayer polycrystalline graphene nanoribbons: effects of interlayer interaction, grain size, and vacancy defects

writerZi-Qiang Yang, Muhammad Asif Shakoori & Hai-Peng Li

Vol.34 (Feb) 2024 | Article no.1 2024


In this paper, molecular dynamics simulations have been employed to investigate the phonon thermal transport in bilayer polycrystalline graphene nanoribbon (pGNR/pGNR), compared with bilayer graphene nanoribbon (GNR/GNR) and pGNR/GNR heterostructure. The interfacial thermal resistance (ITR) of bilayer structures was also calculated using the heat dissipation method. The effects of interlayer interaction, grain size, and vacancy defects on ITR and in-plane phonon thermal conductivity of bilayer structures were investigated. It was found that the ITR as well as in-plane phonon thermal conductivity of pGNR/pGNR was less than that of pGNR/GNR and much less than that of GNR/GNR, for the same size. For the studied bilayer structures, both the ITR and in-plane phonon thermal conductivity decrease with increasing interlayer interactions. Moreover, ITR increases with increasing grain area size whereas decreases with increasing vacancy defects in pGNR-based bilayers. The introduction of pGNR interface roughness and vacancy defects results in an enhanced phonon coupling in pGNR-based bilayers compared to pure GNR/GNR bilayers. Presented simulation investigations will help to understand the interlayer thermal transport properties of polycrystalline graphene and provide essential guidance for experimentally regulating phonon thermal transport between layers of polycrystalline graphene.

1 Introduction

Graphene is a single layer of carbon atoms arranged in a hexagonal lattice structure. It is incredibly thin, lightweight, transparent, strong, excellent conductor of heat and electricity, flexible, and stretchable [1,2,3,4,5]. Graphene-based materials have various shapes and structures, each with unique properties, such as graphene sheets [6], oxide [7], nanoribbons [8], nanotubes [9], and nanoplatelets [10]. Due to its remarkable electronic, mechanical, and thermal properties, it holds immense potential for diverse applications, such as energy storage [2] and conversion [11], electronics and photonics [4], sensors [3], water purification [12], and reinforced composite materials [10]. Numerous properties of graphene-based materials can be enhanced by changing their structural orientations. Thermal transport is one of them, which is very important for improving existing technologies and developing new ones where heat transport and heat management play a crucial role [7, 13,14,15].

Polycrystalline graphene consists of many grains with random orientation, and it can be produced through numerous methods, such as chemical vapor deposition (CVD) on metal substrate, reducing graphene oxide, and annealing process [5, 16]. Graphene nanoribbons (GNRs) have been extensively studied, such as the effects of grain boundaries on charge transport [17], the impact of graphene performance [5, 16], and bilayer graphene has tunable band gaps [18]. Polycrystalline graphene overcomes the disadvantage of single-layer graphene lacking band gaps. Especially in 2018, Cao et al. discovered that superconductivity was observed in twisted bilayer graphene when the twist angle was close to 1.1° [19, 20], which sparked interest in the study of this special structure. The bilayer GNRs are also intriguing due to their potential to blend the electronic properties of GNRs with the unique features of bilayer graphene [21]. Reducing interfacial thermal resistance (ITR) is crucial in applications like thermal interface materials, where efficient heat dissipation is necessary to prevent electronic devices from overheating. ITR is the main obstacle to heat flow from one layer to another in multi-layer materials [1]. With the penetration of semiconductor technology into other applications, the problem of ITR has been encountered in many new fields and devices, such as light-emitting diodes [22], quantum cascade lasers [23], phase change memories [24], thermoelectric devices [11], wearable devices [25], photovoltaic cells [26], and three-dimensional architecture chips [27]. Graphene is a good heat-dissipation material and has ultra-high thermal conductivity. Still, due to its two-dimensional structure limitations, its out-of-plane thermal conductivity is far inferior to in-plane thermal conductivity [28]. Therefore, how to improve its out-of-plane thermal conductivity to reduce the ITR of devices has become a hot topic of concern.

Researchers explore various techniques to minimize ITR, such as introducing intermediate layers with compatible atomic structures, optimizing chemical bonding, and controlling surface roughness [29]. These efforts aim to enhance the overall thermal performance of graphene-based materials in practical applications, making them more effective for thermal management and other heat-related technologies [30]. During the process of rotating bilayer single-crystal graphene, researchers found that the ITR increased with the rotation angle from 0° to 30° [13], and the in-plane thermal conductivity first decreased and then increased from 0° to 30° [31]. The thermal transport of polycrystalline graphene has been extensively investigated through theoretical calculations [32, 33], molecular dynamics (MD) simulations [7, 15, 30, 34,35,36], and CVD [14, 29] from different perspectives. For example, Krasavin et al. [33] predicted the effects of grain boundary structure, size, and shape on the phonon thermal conductivity of polycrystalline graphene. MD simulations were performed to compute the effects of grain size, alignment, and temperature by Wang et al. [34], grain boundary and grain-sized by Liu et al. [13], grain size, tensile strain, and temperature by Wu et al. [15], fatigue by Zhang et al. [36] and recently grain, vacancy and strain by Li et al. [30] on phonon thermal transport of polycrystalline graphene. Mu et al. calculated the thermal conductivity of oxidized polycrystalline graphene [7]. In addition, researchers have attempted to regulate the thermal transport properties of bilayer graphene by doping carbon atoms between layers [37]. Previous literature has reported on the frictional mechanics of multi-layer polycrystalline graphene [38] and on the effect of layer number on the thermal transport of polycrystalline graphene [29]. However, there are few reports on the effects of interlayer forces, grain size, and vacancy defect on the thermal transport of bilayer polycrystalline graphene, and thus the related mechanism is still unclear.

In this paper, the in-plane and out-of-plane thermal transport properties of the polycrystalline graphene nanoribbon (pGNR) bilayer and GNR bilayer are calculated by using the MD simulation method, respectively. The effects of interlayer forces, grain size, and vacancy defects on ITR and in-plane thermal conductivity are discussed for pGNR and GNR bilayers. The impact and mechanism of vacancy defects on the ITR of the pGNR bilayer are addressed from the perspective of phonon transport. This study contributes to a deeper understanding of the physical mechanism of thermal transport in the pGNR bilayer. It provides essential guidance for experimental research on regulating the phonon thermal conductivity and ITR in polycrystalline graphene layered materials.

2 Models and methods

2.1 Bilayer pGNR model and structure relaxation

The structure model of monolayer pGNR is formed by the Thiessen polygon method [30, 39] with the open-source software Atomsk [40]. Upon acquiring a model of monolayer pGNR, adjustments were made to the interlayer spacing between the two layers, setting it precisely at d = 0.34 nm, see Fig. 1a. The MD simulations were implemented with the large-scale atomic/molecular massively parallel simulator (LAMMPS) package [41]. The AIREBO potential [42] was applied to describe the C-C interactions in the pGNR and GNR monolayers, and this potential has widely been used to calculate the thermal conductivity of graphene-based nanomaterials [43, 44]. The interlayer Van der Waals interactions are modeled by the Lennard-Jones potential [38]. The relevant parameters of the Lennard-Jones potential are taken from some commonly used in previous calculations for bilayer graphene [13]. In order to ensure stability and computational efficiency, a time step of 0.5 fs was used in the MD simulation to calculate the required properties. The periodic boundary conditions are applied on the layer plane (x and y directions). The pGNR monolayer/bilayer model underwent a comprehensive relaxation process as follows: initially heating up to 1000K over the duration of 0.5 ns in the NPT ensemble. Subsequently, a cooling phase took place for another 0.5 ns, bringing the system down to 300 K. It was followed by a 1.0-ns relaxation period in both the 300 K NPT and NVT ensembles and, ultimately, an additional 1.0-ns relaxation phase in the NVE ensemble to achieve the desired equilibrium state.

Fig. 1
figure 1

Schematics of bilayer models: a pGNR/pGNR, b pGNR/GNR, and c GNR/GNR

2.2 Calculation method of ITR

The ITR is estimated by using the transient pump–probe method, as described in the reference [45]. This method has been widely used to investigate ITR in diverse two-dimensional heterostructures, as demonstrated in previous studies [46,47,48]. The MD simulation imitates the process of the transient pump–probe experiment, also called the heat dissipation method [45]. The simulation process is as follows: the two-layer heterostructure model undergoes complete relaxation initially at room temperature. After attaining the NVE equilibrium state, a thermal pulse of 50.0 ps is applied to one of the layers. During this process, the temperature of the heating layer (T1) rapidly increases to 400 K while the temperature of the undisturbed layer (T2) remains unchanged at 300 K within the NVT ensemble. Following the heat pulses of 50.0 ps, the thermostat is removed, and heat energy dissipates through interface thermal conduction until the system attains a new thermal equilibrium, as shown in Fig. 2a. Throughout the relaxation process, the temperature evolution of each layer and the total energy of the heat-applied layer are recorded at each time step. Due to the vacuum in the surrounding environment, the only thermal channel between monolayers is the exchange of heat energy between interfaces. The energy dissipation change of the heating layer over time is described as

\(\frac{\partial {E}_t}{\partial t}=A\times \left(\frac{T_1-{T}_2}{R}\right)\)

where Et is the total energy of the heating layer during the time t, R is the ITR, A is the effective heat transfer area between the two layers, and T1/T2 is the temperature of the heating layer/undisturbed layer, respectively. In order to reduce data fluctuations, the collected data points are averaged once every 100 steps. The energy Et after integration becomes as,

\({E}_t={E}_0+\left(\frac{A}{R}\right).{\int}_0^t\left(T{}_1-{T}_2\right) dt\)

where E0 represents the initial energy of the heating layer. As shown in Fig. 2b, the ITR is obtained from the linear fitting relationship of Et and the time integral of the temperature difference (T1-T2). The initial data for 100 ps was excluded from linear fitting due to the significant energy disturbance caused by the thermal pulse in the early stages of the thermal relaxation process.

Fig. 2
figure 2

a The evolution of heating layer temperature (T1), undisturbed layer temperature (T2), and heating layer energy (E1) versus time t during the heat dissipation process. b The relationship between the total energy of the heating layer and the time integral of (T1-T2)

2.3 Phonon density of state analysis

In solid-state physics, phonons are the quantized form of lattice vibration, and the room-temperature thermal conduction of crystals is mainly determined by the interaction of phonons. Phonon density of state (PDOS) is an important parameter that reflects the thermal properties of crystals at the microscopic level. The degree of overlap of PDOS of layered materials determines the phonon transmission behavior of the interface. Based on MD simulations, PDOS can be obtained by Fourier transform of the velocity autocorrelation function as [49]

\(D\left(\omega \right)={\int}_0^t\Gamma (t)\exp \left(-\omega t\right) dt\)

where ω is the frequency, D(ω) is the frequency of PDOS, Г(t) is the velocity autocorrelation function, written as Γ(t) = 〈v(t) ⋅ v(0)〉. The thermal resistance of the interface between heterostructure interfaces is closely related to the overlap factor S of PDOS on both sides of the interface [50]

\(S=\frac{\int_{-\infty}^{+\infty }{D}_1\left(\omega \right){D}_2\left(\omega \right) d\omega}{\int_{-\infty}^{+\infty }{D}_1\left(\omega \right) d\omega {\int}_{-\infty}^{+\infty }{D}_2\left(\omega \right) d\omega},\)

where the D1 and D2 are the PDOS of the materials on both sides of the interface. The larger the S value of the overlap factor, the stronger the phonon coupling between the two layers, and thus the more conducive it is to phonon transport.

3 Results and discussion

3.1 Effect of interlayer interaction

In the stable state of two-dimensional layered materials, interlayer forces, commonly referred to as van der Waals forces, exist between the layers. The Lennard-Jones potential is frequently used in MD simulations to describe the van der Waals interactions between atoms. The van der Waals interaction can be characterized by a potential functionU(r) = χV(r), where V(r) is the Lennard-Jones potential in the initial state and χ is the proportional coefficient used to adjust the strength of the interaction. In order to analyze disparities between pGNR and GNR, three different bilayer structures were considered: pGNR/pGNR, pGNR/GNR, and GNR/GNR, as shown in Fig. 1. This selection allowed for observing how varying interlayer interaction strengths influence in-plane phonon thermal transport and the ITR within the bilayer system. Here the GNR and pGNR models with dimensions of 40 nm × 10 nm are specifically selected, and the average grain area size of pGNR is 50.0 nm2. A uniform χ value within the collapse limit was consistently applied to ensure the system’s stability amid fluctuating interlayer interaction strengths.

Figure 3 illustrates that the ITR of the GNR/GNR bilayer is higher than that of the pGNR/GNR heterostructure under the same χ. Moreover, the ITR of the pGNR/pGNR bilayer is the lowest. For example, at χ = 1, the ITR of pGNR/GNR reduces to 64% of the ITR of GNR/GNR. Similarly, the ITR of pGNR/pGNR reduces to 55% of the ITR of GNR/GNR. This indicates that the grain boundaries of pGNR enhanced phonon transport at the interface, leading to significantly improved out-of-plane thermal conduction in the bilayers. The decrease in ITR can be explained by the increased friction caused by the enhancement of interface roughness and defects, resulting in an enhancement of phonon coupling. When the interlayer interaction strength increases as χ = 2, 4, 6, 8, and 10, the ITR of GNR/GNR bilayer decreases to 90%, 74%, 62%, 56%, and 46%, respectively, of the original value (χ = 1). The ITR of pGNR/GNR heterostructure decreases to 92%, 72%, 67%, 54%, and 53% of the original value (χ = 1), respectively. The ITR of the pGNR/pGNR bilayer is reduced to the original 78%, 55%, 48%, 39%, and 28%, respectively. Obviously, the ITR of the pGNR/pGNR bilayer is more susceptible to interlayer forces.

Fig. 3
figure 3

The variation of ITR (R) with the intensity of interlayer interaction (χ)

Figure 4 displays the influence of interlayer interaction strength on the in-plane phonon thermal conductivity (κ) of studied bilayer structures. At χ = 1, our presented κ value for the GNR/GNR bilayer was compatible with earlier reported work [9, 51]. With the rise in χ, the in-plane κ of the bilayer structures decreases to varying degrees. The in-plane κ of pGNR/pGNR decreased to the original 94%, 91%, 84%, 82%, and 76%, while the in-plane κ of GNR/GNR decreased to the original 95%, 91%, 88%, 84%, and 80%, respectively, for χ = 2, 4, 6, 8, and 10. Interestingly, the in-plane κ of pGNR/GNR is lower than that of GNR/GNR and higher than that of pGNR/pGNR at the same χ. This phenomenon arises due to the distinct grain orientations and defects in various grain domains within pGNR. All in all, introducing interlayer interaction is analogous to imposing strain in the out-of-plane direction, thereby influencing the overall structure of the model. The larger interlayer interaction intensity corresponds to the larger van der Waals interaction which reduces in-plane phonon thermal conductivity because of suppressing the contribution of the out-of-plane acoustic mode (ZA) [52].

Fig. 4
figure 4

The variation of in-plane thermal conductivity (κ) with the intensity of interlayer interaction (χ). Here the GNR and pGNR models with dimensions of 40 nm × 10 nm are specifically selected, and the average grain area size of pGNR is 50.0 nm2

3.2 Effect of grain size

Figure 5 shows the pGNR grain size effect on the ITR of pGNR-based bilayers. Three randomly generated models with the same grain size were selected for simulation, and the average value was taken to reduce computational fluctuations. Upon reaching a dimension of 400 nm2, corresponding to a bilayer GNR model, the ITR measures 2.82 × 10−2 m2K/MW, aligning with findings reported in the relevant literature [13]. When the average grain size of pGNR is 25, 50, 100, and 200 nm2, the ITRs of the pGNR/pGNR bilayer (pGNR/GNR bilayer) decrease to be 55% (60%), 62% (67%), 74% (89%), 91% (96%) of ITR of GNR/GNR, respectively. The surface of GNR is smooth and complete, resulting in lower phonon coupling strength compared to pGNR in bilayer structures. As the grain size of pGNR increases, the corresponding grain boundaries and defects decrease. Thus weakening phonon coupling between two layers may be unfavorable for phonon thermal transport at the interface.

Fig. 5
figure 5

The variation of ITR (R) of bilayer structures with the average grain area size of pGNR layer. Here the GNR and pGNR models with dimensions of 40 nm × 10 nm are specifically selected

3.3 Effect of vacancy defect

Due to the inevitable occurrence of vacancy defects in layer materials in practical applications, it is important to study the influence of vacancy defects on thermal transport at interlayer interfaces. Figure 6 shows the ITRs of the GNR/GNR bilayer and pGNR/pGNR bilayer under the influence of random vacancy defects rate (ε), respectively. Here vacancy defects rate is generally defined as the ratio between the number of vacancies and the number of atoms in a perfect crystal. From Fig. 6, it can be seen that the ITR decreases with the increase of ε both in the GNR/GNR bilayer and in the pGNR/pGNR bilayer. However, the sensitivity of the ITR to ε is inconsistent during these two bilayer models. As the ε increases by 1%, 3%, and 5%, the ITR of GNR/GNR bilayer (pGNR/pGNR bilayer) is 89% (97%), 61% (85%), and 52% (68%), of the initial state (ε = 0), respectively.

Fig. 6
figure 6

The variation of ITR (R) of GNR/GNR and pGNR/pGNR bilayer structures versus random vacancy defect rate (ε). Here the bilayer models with dimensions of 40 nm × 10 nm are specifically selected, and the average grain area size of pGNR is 50.0 nm2

Figure 7 shows the out-of-plane PDOS (PDOSz) of the pGNR/pGNR bilayer under different vacancy defect rates ε. It is found that when the ε is 0, 1%, 3%, and 5%, the corresponding overlap factors S are 0.0036, 0.0041, 0.0051, and 0.0052, respectively. The increase in overlap factor S indicates that the increase in vacancy defect rate enhances the heat transfer at the interface, leading to a decrease in ITR. From a structural perspective, as the number of vacancy defects increases, the horizontal binding of atoms decreases while the out-of-plane vibration increases.

Fig. 7
figure 7

The out-of-plane PDOS (PDOSz) of pGNR/pGNR bilayer under different vacancy defect rates (ε)

4 Conclusions

In this paper, MD simulations have been performed to calculate the in-plane phonon thermal conductivity and ITR of pGNR-based bilayers. The effects of interlayer forces, grain size, and vacancy defects on phonon thermal transport were studied systematically. It was found that the in-plane phonon thermal conductivity of pGNR is lower than that of GNR with the same size, regardless of whether it is a single layer or a bilayer, and weakens with the increase of interlayer forces in the bilayer structure. This is because the interface distance between bilayer membranes decreases with the increase of interlayer interaction intensity, and the contribution of phonon bending mode (ZA) is suppressed. More significant interlayer interaction intensity corresponds to a larger van der Waals interaction, which reduces the in-plane phonon thermal conductivity. The ITR of the bilayer structure will decrease with the increase of interlayer force, the decrease of crystal grain size, and the increase of vacancy defects. The ITR of the pGNR bilayer is smaller than that of the GNR bilayer, which can be explained by the increased friction caused by the enhancement of interface roughness and defects, resulting in an enhanced phonon coupling in the pGNR-based bilayers, compared to the GNR bilayer. The research results of this paper help to understand the interlayer thermal transport properties of pGNR and also provide meaningful guidance for the regulation of phonon thermal transport between layers of pGNR.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.


  1. J. Chen, X. Xu, J. Zhou, B. Li, Interfacial thermal resistance: past, present, and future. Rev. Mod. Phys. 94(2), 025002 (2022). https://doi.org/10.1103/RevModPhys.94.025002

  2. D.A.C. Brownson, D.K. Kampouris, C.E. Banks, An overview of graphene in energy production and storage applications. J. Power Sources 196(11), 4873–4885 (2011). https://doi.org/10.1016/j.jpowsour.2011.02.022

  3. Y. Shao, J. Wang, H. Wu, J. Liu, I.A. Aksay, Y. Lin, Graphene based electrochemical sensors and biosensors: a review. Electroanal. 22(10), 1027–1036 (2009). https://doi.org/10.1002/elan.200900571

  4. P. Avouris, F. Xia, Graphene applications in electronics and photonics. MRS. Bullet. 37(11), 1225–1234 (2012). https://doi.org/10.1557/mrs.2012.206

  5. C.H. Lui, Z. Li, K.F. Mak, E. Cappelluti, T.F. Heinz, Observation of an electrically tunable band gap in trilayer graphene. Nat. Phys. 7(9), 944–947 (2011). https://doi.org/10.1038/nphys2102

  6. B. Wang, Q. Cao, W. Shao, Z. Cui, Effect of vacancy defects on the heat transfer coefficient of partially stacked graphene sheets. J. Mater. Sci. 57(4), 8167–8179 (2022). https://doi.org/10.1007/s10853-022-07128-5

  7. X. Mu, Z. Song, Y. Wang, Z. Xu, D.B. Go, T. Luo, Thermal transport in oxidized polycrystalline graphene. Carbon. 108(11), 318–326 (2016). https://doi.org/10.1016/j.carbon.2016.07.023

  8. Y. Xu, X.Y. Xu, W. Zhang, T. Ouyang, C. Tang, Thermoelectric properties of polycrystalline graphene nanoribbons. Acta. Phys. Sin. 68(24), 247202 (2019). https://doi.org/10.7498/aps.68.20191276

  9. C. Diao, Y. Dong, J. Lin, Reactive force field simulation on thermal conductivities of carbon nanotubes and graphene. Int. J. Heat. Mass. Transf. 112(9), 903–912 (2017). https://doi.org/10.1016/j.ijheatmasstransfer.2017.05.036

  10. D. Chen, G. Chen, In situ synthesis of thermoplastic polyurethane/graphene nanoplatelets conductive composite by ball milling. J. Reinf. Plast. Comp. 32(5), 300–307 (2013). https://doi.org/10.1177/0731684412471230

  11. W. Liu, S. Bai, Thermoelectric interface materials: a perspective to the challenge of thermoelectric power generation module. J. Materiomics. 5(3), 321–336 (2019). https://doi.org/10.1016/j.jmat.2019.04.004

  12. B. Wang, Q. Liu, Z. Fan, A mini review: application progress of magnetic graphene three-dimensional materials for water purification. Front. Chem. 8(11), 595643 (2020). https://doi.org/10.3389/fchem.2020.595643

  13. W. Liu, Y. Hong, J. Zhang, Y. Yue, Anisotropic Thermal Transport in Twisted Bilayer Graphene. Phys. Chem. Chem. Phys. 24(8), 21722 (2022). https://doi.org/10.1039/D2CP01896A

  14. Z. Aksamija, I. Knezevic, Lattice Thermal Transport in Large-Area Polycrystalline Graphene. Phys. Rev. B. 90(7), 035419 (2014). https://doi.org/10.1103/PhysRevB.90.035419

  15. P.H. Wu, S.S. Quek, Z.D. Sha, Z.L. Dong, X.J. Liu, G. Zhang, Q.X. Pei, Y.W. Zhang, Thermal transport behavior of polycrystalline graphene: a molecular dynamics study. J. Appl. Phys. 116(11) (2014). https://doi.org/10.1063/1.4902852

  16. S. Ghosh, W. Bao, D.L. Nika, S. Subrina, E.P. Pokatilov, C.N. Lau, A.A. Balandin, Dimensional crossover of thermal transport in few-layer graphene. Nat. Mater. 9(5), 555–558 (2010). https://doi.org/10.1038/nmat2753

  17. J. Wu, Y. Li, D. Pan, C. Jiang, C. Jin, F. Song, G. Wang, J. Wan, Effect of grain boundaries on charge transport in CVD-grown bilayer graphene. Carbon. 147(6), 434–440 (2019). https://doi.org/10.1016/j.carbon.2019.03.029

  18. Y. Zhang, T.T. Tang, C. Girit, Z. Hao, M.C. Martin, A. Zettl, M.F. Crommie, Y.R. Shen, F. Wang, Direct observation of a widely tunable bandgap in bilayer graphene. Nature. 459(6), 820–823 (2009). https://doi.org/10.1038/nature08105

  19. Y. Cao, V. Fatemi, A. Demir, S. Fang, S.L. Tomarken, J.Y. Luo, J.D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R.C. Ashoori, P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices. Nature. 556(3), 80–84 (2018). https://doi.org/10.1038/nature26154

  20. Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices. Nature. 556(3), 43–50 (2018). https://doi.org/10.1038/nature26160

  21. M.H. Oliveira, J.M.J. Lopes, T. Schumann, L.A. Galves, M. Ramsteiner, K. Berlin, A. Trampert, H. Riechert, Synthesis of quasi-free-standing bilayer graphene nanoribbons on SiC surfaces. Nat. Commun. 6(7), 7632 (2015). https://doi.org/10.1038/ncomms8632

  22. X. Lei, H. Zheng, X. Guo, Z. Zhang, J. Wu, C. Xu, S. Liu, Reduction of die-bonding interface thermal resistance for high-power LEDs through embedding packaging structure. IEEE. Trans. Power Electron. 32, 5520 (2017). https://doi.org/10.1109/ITHERM.2016.7517664

  23. S. Takagi, H. Tanimura, T. Kakuno, R. Hashimoto, K. Kaneko, S. Saito, Evaluation of heat dissipation characteristics of quantum cascade laser with diamond submount using structure function and three-dimensional thermal fluid simulation. Jpn. J. Appl. Phys. 60(12), 124003 (2021). https://doi.org/10.35848/1347-4065/ac33d0

  24. J.P. Reifenberg, K.W. Chang, M.A. Panzer, S. Kim, J.A. Rowlette, M. Asheghi, H.S.P. Wong, K.E. Goodson, Thermal boundary resistance measurements for phase-change memory devices. IEEE. Electron. Device. Lett. 31(1), 56 (2010). https://doi.org/10.1109/LED.2009.2035139

  25. Y.S. Ju, Thermal management and control of wearable devices. iScience. 25(7), 104587 (2022). https://doi.org/10.1016/j.isci.2022.104587

  26. J. Zhang, H. Zhai, Z. Wu, Y. Wang, H. Xie, M. Zhang, Enhanced performance of photovoltaic–thermoelectric coupling devices with thermal interface materials. Energy. Rep. 6(11), 116–122 (2020). https://doi.org/10.1016/j.egyr.2019.12.001

  27. M.M. Waldrop, The semiconductor industry will soon abandon its pursuit of Moore’s law. Nature. 530(2), 7589 (2016). https://doi.org/10.1038/530144a

  28. Z. Fan, L.F.C. Pereira, P. Hirvonen, M.M. Ervasti, K.R. Elder, D. Donadio, T. Ala-Nissila, A. Harju, Thermal conductivity decomposition in two-dimensional materials: application to graphene. Phys. Rev. B. 95(4), 144309 (2017). https://doi.org/10.1103/PhysRevB.95.144309

  29. D. Estrada, Z. Li, G.M. Choi, S.N. Dunham, A. Serov, J. Lee, Y. Meng, F. Lian, N.C. Wang, A. Perez, R.T. Haasch, J.M. Zuo, W.P. King, J.A. Rogers, D.G. Cahill, E. Pop, Thermal transport in layer-by-layer assembled polycrystalline graphene films, NPJ 2D. Mater. Appl. 3(3), 10 (2019). https://doi.org/10.1038/s41699-019-0092-8

  30. Z. Yang, R. Wang, H. Li, H.K. Tang, K. Han, Phonon thermal transport in polycrystalline graphene: effects of grain, vacancy and strain. Int. J. Heat. Mass. Transf. 209(8), 124057 (2023). https://doi.org/10.1016/j.ijheatmasstransfer.2023.124057

  31. S. Han, X. Nie, S. Gu, W. Liu, L. Chen, H. Ying, L. Wang, Z. Cheng, L. Zhao, S. Chen, Twist-angle-dependent thermal conduction in single-crystalline bilayer graphene. Appl. Phys. Lett. 118(5), 193104 (2021). https://doi.org/10.1063/5.0045386

  32. D.V. Kolesnikov, V.A. Osipov, Low-Temperature Thermal Conductivity in Polycrystalline Graphene. EPL. 100(10), 26004 (2012). https://doi.org/10.1209/0295-5075/100/26004

  33. S.E. Krasavin, V.A. Osipov, Impact of grain boundary characteristics on thermal transport in polycrystalline graphene: analytical results. J. Appl. Phys. 125(3), 084301 (2019). https://doi.org/10.1063/1.5082151

  34. Y. Wang, Z. Song, Z. Xu, Characterizing phonon thermal conduction in polycrystalline graphene. J. Mater. Res. 29(2), 362–372 (2014). https://doi.org/10.1557/jmr.2013.380

  35. H.K. Liu, Y. Lin, S.N. Luo, Grain boundary energy and grain size dependences of thermal conductivity of polycrystalline graphene. J. Phys. Chem. C. 118(10), 24797 (2014). https://doi.org/10.1021/jp508035b

  36. S. Zhang, J. Zhang, Fatigue and its effect on the mechanical and thermal transport properties of polycrystalline graphene. J. Mater. Sci. 56(3), 10367–10381 (2021). https://doi.org/10.1007/s10853-021-05953-8

  37. Y. Mao, G. Malcolm Stocks, J. Zhong, First-principles study of the doping effects in bilayer graphene. New. J. Phys. 12(3), 033046 (2010). https://doi.org/10.1088/1367-2630/12/3/033046

  38. A. Kavalur, W.K. Kim, Molecular dynamics study on friction of polycrystalline graphene. Comput. Mater. Sci. 137(9), 346–361 (2017). https://doi.org/10.1016/j.commatsci.2017.06.006

  39. X. Zheng, T. Sun, J. Zhou, R. Zhang, P. Ming, Modeling of polycrystalline material microstructure with 3d grain boundary based on Laguerre–Voronoi Tessellation Mater. 15(6), 1996 (2022) https://doi.org/10.3390/ma15061996

  40. P. Hirel, Atomsk: a tool for manipulating and converting atomic data files. Comput. Phys. Commun. 197(12), 212–219 (2015). https://doi.org/10.1016/j.cpc.2015.07.012

  41. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117(1), 1–19 (1995). https://doi.org/10.1006/jcph.1995.1039

  42. S.J. Stuart, A.B. Tutein, J.A. Harrison, A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 112(4), 6472–6486 (2000). https://doi.org/10.1063/1.481208

  43. A. Fox, U. Ray, T. Li, Thermal conductivity of graphene grain boundaries along arbitrary in-plane directions: a comprehensive molecular dynamics study. J. Appl. Phys. 125(12), 015101 (2019). https://doi.org/10.1063/1.5059561

  44. C. Wang, S. Lu, X. Yu, H. Li, Alkyl group functionalization-induced phonon thermal conductivity attenuation in graphene nanoribbons. Chin. Phys. B. 28, 016501 (2019). https://doi.org/10.1088/1674-1056/28/1/016501

  45. Y. Hong, J. Zhang, X.C. Zeng, Interlayer thermal conductance within a phosphorene and graphene bilayer. Nanoscale. 8(10), 19211–19218 (2016). https://doi.org/10.1039/c6nr07977f

  46. D. Han, X. Wang, W. Ding, Y. Chen, J. Zhang, G. Xin, L. Cheng, Phonon thermal conduction in a graphene-C3N Heterobilayer using molecular dynamics simulations. Nanotechnol. 30(7), 075403 (2019). https://doi.org/10.1088/1361-6528/aaf481

  47. J. Zhou, H. Li, H.K. Tang, L. Shao, K. Han, X. Shen, Phonon thermal transport in silicene/graphene heterobilayer nanostructures: effect of interlayer interactions. ACS Omega. 7(7), 5844–5851 (2022). https://doi.org/10.1021/acsomega.1c05932

  48. Y. Hong, C. Zhu, M. Ju, J. Zhang, X.C. Zeng, Lateral and flexural phonon thermal transport in graphene and stanene bilayers. Phys. Chem. Chem. Phys. 19(2), 6554–6562 (2017). https://doi.org/10.1039/c6cp08276a

  49. K.E. Eshkalak, S. Sadeghzadeh, F. Molaei, Interfacial thermal resistance mechanism for the polyaniline (C3N)–graphene heterostructure. J. Phys. Chem. C. 124(5), 14316–14326 (2020). https://doi.org/10.1021/acs.jpcc.0c02051

  50. C. Si, X.D. Wang, Z. Fan, Z.H. Feng, B.Y. Cao, Impacts of potential models on calculating the thermal conductivity of graphene using nonequilibrium molecular dynamics simulations. Int. J. Heat. Mass. Transf. 107(4), 450–460 (2017). https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.065

  51. 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, L. Shi, Two-dimensional phonon transport in supported graphene. Science. 328(4), 213–216 (2010). https://doi.org/10.1126/science.1184014

  52. M.Li, M.A. Shakoori, R. Wang, H. Li, Phonon thermal transport at the interfaces of graphene/vertically aligned carbon nanotubes/hexagonal boron nitride sandwich heterostructure. Chinese Phys. Lett. 41(1), 016302 (2024). https://doi.org/10.1088/0256-307X/41/1/016302


This work is supported by the Fundamental Research Funds for the Central Universities of China (Grant No. 2019ZDPY16). We are also grateful for the support of the funding for the key discipline of physics at the China University of Mining and Technology.

Author information

Authors and Affiliations


Conceptualization: Z.Q. and H.P.; methodology and investigation: Z.Q. and M.A.; writing—original draft preparation: Z.Q.; writing—review and editing: M.A. and H.P. All authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Muhammad Asif Shakoori or Hai-Peng Li.

Ethics declarations

Ethics approval and consent to participate

The authors declare they have upheld the integrity of the scientific record.

Consent for publication

The authors give their consent for publication of this article.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

[Source: https://link.springer.com/article/10.1007/s43673-023-00111-5]