ES Materials & Manufacturing, 2019, 3, 16-21

Published online: 29 Jan 2018

Received 12 Dec 2018, Accepted 29 Jan 2018

Phonon Thermal Transport Properties of Graphene Periodically Embedded with Four- and Eight-membered Rings: a Molecular Dynamics Study

Chuanjiang Tang,1 Xiaoxiang Yu,2,3 Gen Li,1 Nuo Yang2,3 and Jingtao Lü1, *

1School of Physics and Wuhan National High Magnetic Field Center, Huazhong University of Science and Technology, 430074 Wuhan, China

2State Key Laboratory of Coal Combustion, Huazhong University of Science and Technology, 430074 Wuhan, China

3Nano Interface Center for Energy (NICE), School of Energy and Power Engineering, Huazhong University of Science and Technology, 430074 Wuhan, China

*E-mail: [email protected] (N. Yang); [email protected] (J. Lu)


It is studied the thermal conductivity of graphene periodically embedded with four- and eight- membered rings (GFERs) by using non-equilibrium molecular dynamics simulations. This kind of structure has been experimentally synthesized recently. The dependence of thermal conductivity on the length (L) and temperature(T) is investigated. It is found that the thermal conductivity of GFERs is significantly lower than that of pristine graphene. On the other hand, the length dependence of thermal conductivity follows ~ logL behavior. In addition, the temperature dependence of thermal conductivity of GFERs follows ~T−α behavior. It is also found that there exists a large thermal rectification (TR) in graphene-GFERs heterostructures, the heat flux from the pristine graphene to the GFERs direction is larger than that in the opposite direction. The dependence of the TR ratio on system parameters is investigated.

Table of Content

Brand new thermal transport properties of graphene periodically embedded with four- and eight membered rings


Keywords: Thermal conductivity     Thermal rectification      Molecular dynamics simulation     Temperature dependence     Size effect     Four and eight-membered rings

1. Introduction

Graphene is made of single atomic layer of carbon atoms and forms a two-dimensional honeycomb lattice.1 As the thinnest two-dimensional material discovered so far, it has excellent thermal, electrical properties, thus attracted enormous interest, since it firstly became experimentally accessible.2 Various numerical and theoretical methods have been used to study the lattice thermal conductivity of different graphene structures, among which are the Boltzmann transport equation (BTE) with the relaxation time approximation,3 Monte Carlo (MC) simulations,4 and Molecular Dynamics (MD) simulations.5–7 MD simulation is a powerful tool to study the thermal properties of nanostructures. Using this approach, thermal conductivity of graphene nanoribbons,8–10 graphene phononic crystal,11–13 graphene disk,14 graphene with impurities and defects15,16 have been studied. Comparing with pristine graphene, the thermal conductivity of these structures is greatly reduced.

Fig. 1 (a) Schematic picture of GFERs. (b) Schematic picture of graphene-GFERs heterostructure. Periodic boundary condition is used along transverse direction and fixed boundary condition is used in longitudinal direction. Heat bath with higher temperature (red box) TL and lower temperature (blue box) TR are applied at two ends.

Embedding non-hexagonal rings into sp2-hybridized carbon networks is considered a promising strategy to enrich the family of low-dimensional graphene structures. However, non-hexagonal rings are energetically unstable compared to the hexagonal counterparts. Very recently, Liu et al.17 reported on-surface synthesis of graphene-like nanoribbons with periodically embedded four- and eight-membered rings. The thermal transport of this kind of structure has not been considered yet.

Recently, the thermal rectification effect attracts more attentions due to its vital role in thermal devices18. The thermal rectification means that the magnitude of heat flux changes when the temperature gradient is reversed in direction. However, it is not easy to obtain a high-performance thermal rectifier. The thermal rectification effect has been found, by theoretical calculations or experimental measurements, in nonlinear lattice models,19 graded mass density,20 asymmetric structure,9,13 quantum system,21,22 interfacial structures,23 or phase change materials24,25. Low dimensional structures tend to realize thermal rectification effect due to ballistic transport of phonons, such as in graphene ribbons9 and graphene phononic crystals.13

In this paper, using non-equilibrium molecular dynamics (NEMD) simulations, the size and temperature dependence of GFERs thermal conductivity have been studied. It is found a drastic drop in the thermal conductivity compared to similar pristine graphene structure. It is also investigated that the thermal rectification property of the graphene-GFERs heterostructure, where the heat flux in one direction is different from that in the opposite direction, and large thermal rectification (TR) ratio is found.

2. Molecular dynamics simulation methods

Fig. 1(a) is a schematic picture of GFERs, and Fig. 1(b) is graphene-GFERs heterostructure where the left region is pristine graphene, and the right region is GFERs. In our calculations, the lattice constant, a, and thickness of graphene are chosen as 1.418 Å and 3.34 Å, respectively. The adaptive intermolecular reactive empirical bond order (AIREBO) potential is used to describe the carbon-carbon interaction, which includes both two-body and three-body potential terms.26–28 The fixed (periodic) boundary condition is adopted along the longitudinal (transversal) direction.29 Non-equilibrium MD simulation is performed using LAMMPS30 with the Langevin approach. The equation of motion for the i-th degrees of freedom that couples to the Langevin bath is written as

                                                                     pidt = - Hqi + ξi - 1mλp        (1)

where H is the Hamiltonian of the system, pi and qi are the momentum and coordinate of particle i, λ is the dissipation rate, and ξi is the Gaussian random force with zero mean and variance 2mλkBT according to the fluctuation-dissipation theorem. There is a free parameter, λ, which controls the coupling between the system under study and heat bath. For those degrees of freedom that do not couple directly to the baths, the equations of motion do not have the last two terms at the right side of the Eq. (1).31,32

Fig. 2 Impacts of heat bath parameter on thermal properties of GFERs. (a) Temperature profile with parameter λ of Langevin heat bath. (b) Heat flux per cross section area versus λ in logarithmic scale for Langevin heat bath. Here the λ is from 0.01 to 50.0.

The Langevin equation is integrated by using the velocity Verlet method with a time step of 0.25 fs. The structures are first relaxed at constant temperature, T , for 3×108 time steps using the Nosé-Hoover thermostat. In order to establish a temperature gradient, the atoms at the two ends are coupled with a Langevin heat bath, whose temperature is TL = T (1+∆) in the hot bath, and TR = T (1−∆) in cold bath.  T is average temperature, 2∆ is the relative temperature bias, and steady state is achieved after 1 ns. The simulation is then continued for another 2 ns for data collection.7 The thermal conductivity κ is calculated based on Fourier’s Law as33,34

                                                                                  κ = - JWd(T/x)      (2)

where J is the heat flux, which is defined as J = (dEhot/dt+dEcold/dt)/2, Ehot and Ecold are the total energy subtracted from (added by) the hot bath (cold bath), respectively. W is the width of the GFERs, d is the thickness of the sample as 3.34 Å, ∂T/∂x is the temperature gradient. In order to quantitatively describe the thermal rectification, the TR ratio is defined as

                           η = JLR -J RL J RL × 100 %                                            (3)

where the JLR is the heat flux along the pristine graphene to the GFERs direction, corresponding ∆ > 0, JRL is the heat flux along the GFERs to the pristine graphene direction, corresponding ∆ < 0.

3.  Results and Discussions

Firstly, it is tested that the effect of λ on the calculated κ. To make little effect on the intrinsic properties of the GFERs, λ can not be too large. At the same time, λ can not be too small either. Or else, the time needed to equilibrate the system would be too long, and there will occur a large temperature jump between heat bath and system. Fig. 2(a) shows the results with different λ between 0.01 and 50.0, with other parameters fixed T  = 300 K, L =11.67 nm, W =5.2 nm. It can see that, if λ is too large, the degrees of freedom that couple to the thermal bath can not equilibrate with the baths. As a results, the temperature gradient in the system is small. Consequently, the heat current JL decreases in Fig. 2(b). This may influence the accuracy of the numerical results. Based on these considerations, an intermediate value of λ = 0.05 is used for Langevin heat baths.

For the Langevin MD, at two ends of the structure, fixed boundary condition is imposed on the boundary atoms. The heat baths are attached to neighboring atoms.35 But the number of atoms that couple directly to the Langevin bath is a parameter to be chosen. Fig. 3 shows the temperature profile of the structure with different number of atomic layers (NL) attached to the Langevin heat bath. Fig. 4(a) and Fig. 4(b) show the dependence of heat flux per cross section area and the calculated thermal conductivity on NL, respectively. It can be seen that the result converges at NL = 4. This is the parameter used in the following calculations.36,37

Fig. 3  Temperature profile with different number of layers where the Langevin heat bath is applied.

Fig. 4  Impacts of number of heat bath layer on thermal properties of GFERs. (a) Heat flux JL per cross section area along the longitudinal direction versus NL. (b) Thermal conductivity κ versus NL.

When the width of simulation cell is not sufficiently large, there will occur a size effect on the thermal conductivity.38 Here, it is studied the dependence of thermal conductivity on the width of the simulation cell. W is changed from 2.6 nm to 7.8 nm. In Fig. 5(a), the results show that with the increase of width, there is little difference between the temperature profile, and the fluctuation gradually becomes smaller when the width is larger.39 The calculated thermal conductivity converges when the width is larger than 5.2 nm in Fig. 5(b). This tendency is consistent with the other study on graphene structure like GNRs,40 GPnC,12 or other graphene structures with defects.7 Hence, in the following simulation, GFERs with 5.2 nm in wide are used to eliminate size effect in the transverse direction.41

Fig. 5  Impacts of width on thermal properties of GFERs. (a) Temperature profile versus width. (b) κ versus width W from 2.6 nm to 7.8 nm. The other parameter L =17 nm,  T=300 K.

After determining the parameters, it is considered that the impact of system length L on the thermal conductivity. Previous study using MD simulation shows that the thermal conductivity increases with L as κ logL.42–44 Here, L is changed from 12 nm to 54 nm at different temperatures. The κ increases from 42.1 W/mK to 92.5 W/mK at 300 K. Results at other temperatures (300 K to 1000 K) are all shown in Fig. 6(a). It exhibits the same behavior as pristine graphene, which has been found in experimental study.42 It can be seen that, κ of GFERs is much smaller than that of pristine graphene. It has been verified that the thermal conductivity contributed by short-range acoustic and optical phonons is temperature independent in previous works.45,46 In graphene, the long-range acoustic phonons dominate the thermal transport.47 The significant decrease of thermal conductivity of GFERs is probably be caused by the affected long-range acoustic phonons due to change of the dispersion relation.29

Fig. 6  (a) Impacts of length on thermal properties of GFERs and pristine graphehe. The length L increase from 12 nm to 54 nm. The parameters are W = 5.2 nm,  T = 300 K for pristine graphene and increases from 300 K to 600 K for GFERs. The thermal conductivity κ versus the length L in log scale. (b) The heat flux versus the length of the graphene-GFERs heterostructure system. The thermal rectification ratio versus the length is shown in insets.

Fig. 6(b) shows the heat flux from two different directions versus the length in the graphene-GFERs heterostructure system, the red line represents the heat flux along the pristine graphene to the GFERs direction, the black is the heat flux from opposite direction, and the insets is the TR ratio versus the system length.  The average temperature T  is fixed as 300 K, the normalized temperature bias ∆ = 0.3, and change the length L from 12 nm to 60 nm, the result show that the thermal rectification decreases with the increasing of the system length, which due to the thermal transport transits to the diffusive regime.7,15 But when the length is more than 60 nm, the TR ratio still have 65%. This indicates that such a heterostructure is a very excellent thermal rectifier which can be suitable for a wide range of system length.

Then, it is discussed the temperature dependence of κ where the length is changed from 12 nm to 27 nm in Fig. 7(a). It is found that κ decrease with the average temperature increasing as Tα12,48, due to the Umklapp phonon-phonon scattering, graphene has a T−1 behavior which is consistent with previous work.49,50 However, for GFERs, the values of α are 0.59, 0.75, 0.76 and 0.78. Obviously, the power exponents (α) increases with increasing the length.51 Thus, it is expected that  α approaches 1 with increase of the length.

Fig. 7  (a) Thermal conductivity of GFERs and pristine graphene with different average temperature T. The T increases from 300 K to 1000 K. The parameters are W=5.2 nm, the length L is 17 nm for grapnene and increase from 12 nm to 27 nm for GFERs. The figure is plotted in a log-log scale. The symbols are numerical data and the lines are fitted lines. (b) The heat flux versus the average temperature T of the graphene-GFERs heterostructure system, and the thermal rectification ratio versus T is shown in insets.

Fig. 8 The dependence of the normalized temperature bias ∆ of the graphene-GFERs heterostructure system on the thermal rectification. (a)The heat flux versus normalized temperature bias ∆, the red line represents the heat flux along the pristine graphene to the GFERs direction, and the black line represents the heat flux from the opposite direction. (b) The thermal rectification ratio versus ∆.

It is also discussed that the temperature dependence of the TR ratio in the graphene-GFERs heterostructure system. L is 48 nm, ∆ is 0.2, and the average temperature T changes from 300 K to 800 K. From Fig. 7(b) it can be found that with the T increasing, the heat flux along two directions both increases, but the growth rate of the heat flux along the pristine graphene to the GFERs direction is smaller than that in the opposite direction. As a result, the calculated TR ratio decreases with the T increasing. The largest TR ratio is 30% in 300 K, but only 1% in 800 K. It is because that thermal rectification is mainly caused by the mismatch of the phonon spectrum on both sides of the rectifier, but increasing the average temperature can weaken this mismatch effect, hence the TR ratio is restrained.5,9

In the following, a discussion on the mechanism of rectification in graphene-GFERs heterostructure system is made based on the difference between thermal resistance of two directions, (TL TR)/JLR and (TR TL)/JRL, which has a reciprocal relationship with the thermal conductivity. As shown in Fig. 7(a), the resistances of GFERs are more than three times larger than those of graphene. So, the resistances of GFERs dominates in the total resistance of heterostructure. The resistances depend on the temperature. The difference between two directions total resistance originate mainly from the resistance of GFERs at different temperatures. When the heat transfers from GFERs to graphene, the GFERs region has a higher temperature, corresponding to a lower heat flux and higher resistance. For the inverse direction, from graphene to GFERs, the GFERs region has a lower temperature, corresponding to a higher heat flux and lower resistance. It means the heat flux preferring from GFERs to graphene. This is different from the situation in graphene-graphene phononic crystal system,13 because the resistance of graphene phononic crystal is independent on the temperature and the resistance of graphene dominates in the total resistance of graphene-graphene phononic crystal system.

Finally, the dependence of the normalized temperature bias ∆ of the graphene-GFERs heterostructure system on the thermal rectification is studied. The average temperature T is fixed as 300 K, the system length L is 48 nm, and the ∆ changes from 0.1 to 0.6. Fig. 8(a) is the heat flux along two directions versus the ∆, it can be seen that with the increasing of the ∆, the heat flux along the pristine graphene to the GFERs direction increases steadily, while the heat flux from the opposite direction increases firstly but then becomes stable. Thus the TR ratio increases with the increasing of the ∆19,52. From Figure 8(b) it can be seen that, when the ∆=0.6, the TR ratio is up to 100%, which exceeds many other graphene based thermal rectifiers.5,7 It indicates that this asymmetric structure is a good thermal rectifier which have a great potential in the field of thermal management.

4. Conclusions

In conclusion, it has been calculated the lattice thermal conductivity of graphene periodically embedded with four- and eight-membered rings as realized in recent experiment, and the thermal rectification of the graphene-GFERs asymmetric structure. It is found that, the length and temperature dependence of the thermal conductivity follow similar trends as the pristine graphene structure. But the magnitude of the thermal conductivity drops by several factors. This low thermal conductivity could be useful in thermoelectric and other applications, where reducing the thermal conductivity may help to increase the performance. It is also found that the heat flux along the pristine graphene to the GFERs direction is larger than that in opposite direction, and the smaller system length, lower average temperature, and higher temperature bias is the optimal condition to obtain the higher TR ratio.


This work was supported by the National Natural Science Foundation of China (Grant No. 21873033, 51576067, 51711540031), National Science Foundation of Hubei Province (2017CFA046). The authors thank the National Supercomputing Center in Shanghai, Tianjin (NSCC-TJ) and China Scientific Computing Grid (ScGrid) for providing assistance in computations.



1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. Firsov, Nature, 2005, 438, 197-200. CrossRef    View Record in Scopus

2. H. Wang, S. Hu, K. Takahashi, X. Zhang, H. Takamatsu and J. Chen, Nat. Commun.,2017, 8, 15843. CrossRef    View Record in Scopus

3. 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. CrossRef    View Record in Scopus

4. S. Chen and D. C. Chrzan, Phys. Rev. B, 2011, 84. CrossRef    View Record in Scopus

5. J. Hu, X. Ruan, Z. Jiang and Y. Chen, Nano lett., 2009, 9, 2730-2735. CrossRef    View Record in Scopus

6. A. Rajabpour, S. M. Vaez Allaei and F. Kowsary, Appl. Phys. Lett., 2011, 99, 051917. CrossRef    View Record in Scopus

7. Y. Wang, S. Chen and X. Ruan, Appl. Phys. Lett., 2012, 100, 163101. CrossRef    View Record in Scopus

8. A. Khan, I. Navid, M. Noshin, H. M. A. Uddin, F. Hossain and S. Subrina, Electronics, 2015, 4, 1109-1124. CrossRef    View Record in Scopus

9. N. Yang, G. Zhang and B. Li, Appl. Phys. Lett. 2009, 95, 033107. CrossRef    View Record in Scopus

10. Z. Guo, D. Zhang and X. G. Gong, Appl. Phys. Lett., 2009, 95, 163103. CrossRef    View Record in Scopus

11. M. Yarifard, J. Davoodi and H. Rafii-Tabar, Comput. Mater. Sci., 2016, 111, 247-251. CrossRef    View Record in Scopus

12. S. Hu, M. An, N. Yang and B. Li, Nanotechnology, 2016, 27, 265702. CrossRef    View Record in Scopus

13. S. Hu, M. An, N. Yang and B. Li, Small, 2017, 13, 1602726. CrossRef    View Record in Scopus

14. X. W. N. Y. Dengke Ma, Hongru Ding and X. Zhang, Int. J. Heat Mass Transfer, 2017, 108, 940. CrossRef    View Record in Scopus

15. K. Yuan, M. Sun, Z. Wang and D. Tang, Int. J. Therm. Sci., 2015, 98, 24-31. CrossRef    View Record in Scopus

16. A. A. Balandin and D. L. Nika, Mater. Today, 2012, 15, 266-275. CrossRef    View Record in Scopus

17. M. Liu, M. Liu, L. She, Z. Zha, J. Pan, S. Li, T. Li, Y. He, Z. Cai and J. Wang, Nat. Commun., 2017, 8, 14924. CrossRef    View Record in Scopus

18. N. Yang, X. Xu, G. Zhang and B. Li, AIP Adv., 2012, 2, 041410. CrossRef    View Record in Scopus

19. J. Lan and B. Li, Phys. Rev. B, 2006, 74, 214305. CrossRef    View Record in Scopus

20. N. Yang, N. Li, L. Wang and B. Li, Phys. Rev. B, 2007, 76, 020301. CrossRef    View Record in Scopus

21. L. Zhang, Y. Yan, C. Q. Wu, J. S. Wang and B. Li, Phys. Rev. B, 2009, 80, 172301. CrossRef    View Record in Scopus

22. L. Zhang, J. S. Wang and B. Li, Phys. Rev. B, 2010, 81, 100301. CrossRef    View Record in Scopus

23. W. Kobayashi, Y. Teraoka and I. Terasaki, Appl. Phys. Lett., 2009, 95, 171905. CrossRef    View Record in Scopus

24. F. A. Herrera, T. Luo and D. B. Go, J. Heat Transfer, 2017, 139, 091301. CrossRef    View Record in Scopus

25. T. Zhang and T. Luo, Small, 2015, 11, 4657-4665. CrossRef    View Record in Scopus

26. D. W. Brenner, O. A. Shenderova, J. Harrison, S. J. Stuart, B. Ni and S. B. Sinnott, J. Phys.: Condens. Matter, 2002, 14, 783-802. CrossRef    View Record in Scopus

27. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472-6486. CrossRef    View Record in Scopus

28. C. Si, X. D. Wang, Z. Fan, Z. H. Feng and B. Y. Cao, Int. J. Heat Mass Transfer, 2017, 107, 450-460. CrossRef    View Record in Scopus

29. L. Yang, J. Chen, N. Yang and B. Li, Int. J. Heat Mass Transfer, 2015, 91, 428-432. CrossRef    View Record in Scopus

30. S. Plimpton, J. Comput. Phys., 1995, 117, 1-19. CrossRef    View Record in Scopus

31. A. Dhar, Adv. Phys., 2008, 57, 457-537. CrossRef    View Record in Scopus

32. J. Chen, G. Zhang and B. Li, J. Phys. Soc. Jpn., 2010, 79, 074604. CrossRef    View Record in Scopus

33. R. Inoue, H. Tanaka and K. Nakanishi, J. Chem. Phys., 1996, 104, 9569-9577. CrossRef    View Record in Scopus

34. X. Yang, A. C. To and R. Tian, Nanotechnology, 2010, 21, 155704. CrossRef    View Record in Scopus

35. G. Hu and B. Cao, Mol. Simul., 2012, 38, 823-829. CrossRef    View Record in Scopus

36. J. W. Jiang, J. Chen, J. S. Wang and B. Li, Phys. Rev. B, 2009, 80, 404. CrossRef    View Record in Scopus

37. R. N. Salaway and L. V. Zhigilei, Int. J. Heat Mass Transfer, 2014, 70, 954-964. CrossRef    View Record in Scopus

38. S. G. Volz and G. Chen, Appl. Phys. Lett., 1999, 75, 2056-2058. CrossRef    View Record in Scopus

39. W. J. Evans, L. Hu and P. Keblinski, Appl. Phys. Lett., 2010, 96, 203112. CrossRef    View Record in Scopus

40. D. Yang, F. Ma, Y. Sun, T. Hu and K. Xu, Appl. Surf. Sci., 2012, 258, 9926-9931. CrossRef    View Record in Scopus

41. H. Wang, K. Kurata, T. Fukunaga, X. Zhang and H. Takamatsu, Int. J. Heat Mass Transfer, 2017, 105, 76-80. CrossRef    View Record in Scopus

42. X. Xu, L. F. Pereira, Y. Wang, J. Wu, K. Zhang, X. Zhao, S. Bae, B. C. Tinh, R. Xie and J. T. Thong, Nat. Commun., 2014, 5, 3689. CrossRef    View Record in Scopus

43. J. Chen, G. Zhang and B. Li, Nanoscale, 2013, 5, 532-536. CrossRef    View Record in Scopus

44. A. Cao, J. Appl. Phys., 2012, 111, 083528. CrossRef    View Record in Scopus

45. A. McGaughey and M. Kaviany, Int. J. Heat Mass Transfer, 2004, 47, 1799-1816. CrossRef    View Record in Scopus

46. A. McGaughey and M. Kaviany, Int. J. Heat Mass Transfer, 2004, 47, 1783-1798. CrossRef    View Record in Scopus

47. T. Feng, X. Ruan, Z. Ye and B. Cao, Phys. Rev. B, 2015, 91, 224301. CrossRef    View Record in Scopus

48. Z. Wang, R. Xie, C. T. Bui, D. Liu, X. Ni, B. Li, and J. T. Thong, Nano Lett., 2011, 11, 113-118. CrossRef    View Record in Scopus

49. J. U. Lee, D. Yoon, H. Kim, W. L. Sang and H. Cheong, Phys. Rev. B, 2011, 83, 081419. CrossRef    View Record in Scopus

50. Y. Kuang, L. Lindsay, S. Shi, X. Wang and B. Huang, Int. J. Heat Mass Transfer, 2016, 101, 772-778. CrossRef    View Record in Scopus

51. N. Yang, S. Hu, D. Ma, T. Lu and B. Li, Sci. Rep., 2015, 5, 14878. CrossRef    View Record in Scopus

52. B. Hu and L. Yang, Chaos: Interdiscip. J. Nonlinear Sci., 2005, 15, 015119. CrossRef    View Record in Scopus