In-plane thermal conductivity of graphene nanomesh: A molecular dynamics study

Computational Materials Science 111 (2016) 247–251

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

In-plane thermal conductivity of graphene nanomesh: A molecular dynamics study M. Yarifard a, J. Davoodi a, H. Rafii-Tabar b,c,⇑ a

Department of Physics, University of Zanjan, Zanjan 45195-313, Iran Department of Medical Physics and Biomedical Engineering, Faculty of Medicine, Shahid Beheshti University of Medical Sciences, Tehran, Iran c Computational Physical Sciences Research Laboratory, School of Nano-Science, Institute for Research in Fundamental Sciences (IPM), Tehran 19395-5531, Iran b

a r t i c l e

i n f o

Article history: Received 30 June 2015 Received in revised form 8 September 2015 Accepted 12 September 2015 Available online 5 October 2015 Keywords: Molecular dynamics Thermal conductivity Graphene nanomesh

a b s t r a c t We have computed the thermal conductivity of the so-called graphene nanomeshes (GNMs) via the Green–Kubo method, employing molecular dynamics (MD) simulation method within a constant(NVE) ensemble. Our results indicate that the thermal conductivity of GNMs is significantly reduced compared with that of pristine graphene so that GNMs may be promising materials in thermoelectric applications. Our simulations exhibit a decreasing behavior and the final convergence of the thermal conductivity as a function of periodicity for a given neck width. Furthermore, our results also show that the thermal conductivity of GNM decreases (increases) as a function of porosity (periodicity) for fixed periodicity (porosity). The effect of the neck width on the thermal conductivity is more pronounced than that of other parameters. It seems that phonon trapping occurs and that the group velocity of heat-carrying phonons decreases. The influence of the geometry of nanoholes on the GNM was also investigated and it was found that GNMs with triangle and square nanoholes, having the same area, offer the lowest and highest thermal conductivities, respectively. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Understanding the physical properties of graphene has been the subject of intense interest and has led to many advances in the fields of nanoelectronics and thermal management ever since the pioneering works of Novoselov et al. [1]. Many potential engineering applications are foreseen for graphene in the realm of thermal management because of its excellent thermal conductivity properties [2–4]. Different forms of graphene have been studied for their novel and diverse applications in transistors [5], sensors [6,7], solar cells [8], antibacterial materials [9] and for energy storage [10]. Both ultra-high thermal conductivity for heat sinking applications [11–13] and ultra-low thermal conductivity for thermoelectric applications [14–16] are associated with this material. Due to strong in-plane sp2 bonding and the low mass of carbon atoms, graphene has a high thermal conductivity along the in-plane directions [17]. Thus, graphene can be one of the best candidates for solving heat dissipation problems in micro/nano electronic devices [2,3]. In principle, the thermal conductivity and the electrical ⇑ Corresponding author at: Computational Physical Sciences Laboratory, School of Nano-Science, Institute for Research in Fundamental Sciences (IPM), Tehran 193955531, Iran. Tel.: +98 2122835058. E-mail address: [email protected] (H. Rafii-Tabar). http://dx.doi.org/10.1016/j.commatsci.2015.09.033 0927-0256/Ó 2015 Elsevier B.V. All rights reserved.

conductivity may be independent in semiconducting nanostructures since different length scales are involved with phonon and electric charges. Electrical conductivity is less sensitive to a decrease in the nanostructure size [16]. Therefore, controlling materials independently of their electrical conductivity is a goal for researchers working on thermoelectric materials for use in energy applications [14]. Defects in the lattice, such as vacancies or Stone–Wales defects [18–20], chemical functionalization [21], strain [22], isotopic impurities (13C) [23], substitutional defects [24] and edge roughness [13,25] can all lead to a reduction in the thermal conductivity of graphene by an order of magnitude below its intrinsic value. Also, grain boundaries in graphene cause the appearance of Kapitza resistance and thermal rectification [26,27]. Therefore, it is plausible that, purposeful engineering of defects can be utilized to accurately tune the thermal properties of graphene. Very recently, following the seminal work of Bai and co-workers on a new graphene nanostructure, referred to as graphene nanomeshes (GNMs) [28], it has been shown that this material has enormous potentials for various applications. Basically, a graphene nanomesh is a piece of graphene sheet containing a periodic array of nanoscale holes. The centre-to-centre distance between two neighbouring nanoholes is called the periodicity, and the smallest edge-to-edge distance between them

248

M. Yarifard et al. / Computational Materials Science 111 (2016) 247–251

is called the neck width, and these form the control parameters for the overall electronic properties of the GNM. Graphene nanomeshes have attracted a great deal of attention, owing to their wide ranging applications in the field of electronics and photonic devices, energy storage, sensors and anti-bacterial materials [10,29–31]. Akhavan and Ghaderi [32] have reported that the graphene nanomesh promises an extremely efficient material for in vivo photothermal therapy in treating cancerous tumours. GNM can be designed with different pore sizes and shapes to filter different gas molecules [33,34] and can also be used for hydrogen storage [35]. Hu [36] has reported that the thermal conductivity of GNMs is an exponential function of the neck width and is independent of the pore lattice arrangement and the pore edge passivation with H atoms for a given density and pore morphology. It is shown that the elastic modulus of GNMs with circular pores scales with the square of density [37]. Also, the thermal expansion of holes in GNMs contradicts the classical prediction [38]. GNM containing a specific hole can tune between metallic and semiconducting states, as well as a semimetal state [39]. Notwithstanding the illuminating insight obtained so far into the electronic properties of GNMs, their thermal properties, namely their thermal conductivities are, however, poorly understood [36]. Since at room temperature the contribution of electrons to the thermal conductivity of graphene is less than that of phonons, classical Molecular Dynamics (MD) simulation [40,41] has been applied to compute the thermal conductivity of graphene as phonons are associated with the displacement of atoms [20,42]. The goal of the present paper is to compute the thermal conductivity of GNMs and its dependence on the periodicity, porosity, neck width and the geometry of the holes. The paper is organized as follows. In Section 2, the computational methodology and the simulation details for calculating the thermal conductivity of GNM are provided. Section 3 reports the simulation results. Finally, concluding remarks of our study are summarized in Section 4.

in-plane heat flux vector ~ J and the heat flux autocorrelation func~ ~ tion hJðtÞ  Jð0Þi for 3–10 ns, depending on the system size to obtain a converged thermal conductivity value. The heat flux vector was computed as [43]

1 ~ JðtÞ ¼

( X

ei~ vi þ

X

i

X 1X ~ ~ v i Þ þ ~rij ð~F j ðijkÞ  ~ v jÞ r ij ðF ij  ~ 2 i;j;i–j i;j;k

)

ð1Þ

where X is the system’s volume calculated as the product of graphene’s area and its thickness that is typically assumed to be equal v and e denote to the graphite interlayer spacing, i.e., 3:35 Å [48], ~ the velocity and energy of atom i, and ~ r and ~ F are the distance and the two/three-body interactions between the carbon atoms. The heat flux was recorded every 5 time steps (2.5 fs) to obtain the heat flux autocorrelation function. The thermal conductivity was obtained via the Green–Kubo theory,

jxx ¼

Z

X kB T

2

Tm

hJxðtÞJxð0Þidt

ð2Þ

0

where kB is the Boltzmann constant and T is the temperature of the system. The bracket term denotes the heat flux autocorrelation function in the x-direction and the average is taken over the time origins. The integration is from zero to a cut-off time (Tm) which is determined by ‘‘first avalanche method” to discard the effect of noise [49]. The final result is the mean value of twenty realizations with different initial conditions. Analogously, we can obtain the thermal conductivity in the y-direction. All of GNMs that we simulated were composed of regular periodic arrays of nanoholes, and the porosity was defined as the ratio of the volume of a nanohole to the volume of the supercell. The geometry of a typical GNM is shown in Fig. 1. We investigated the influence of the nanoholes with different geometries, having the same area, as shown schematically in Fig. 2.

2. Computational methodology and simulation details

3. Simulation results

There are two MD-based approaches for the computation of thermal conductance, namely the equilibrium MD (EMD) and the non-equilibrium MD (NEMD) [43] methods. Thermal conductivities computed via NEMD are commonly smaller than the experimental results, but follow the same trend. EMD simulations have much smaller size-effect in contrast to NEMD, in which heat sources and sinks are employed that scatter phonons. Therefore, EMD is employed in this work. All simulations were carried out using the LAMMPS software package [44,45]. A previous study based on the Boltzmann transport equation has shown that the use of original Tersoff and Brenner potentials underestimate the thermal conductivity of graphene due to the inaccurate handling of phonon dispersion [42]. Therefore, the optimized Tersoff potential is employed in this paper to model the energetics and dynamics of graphene atoms. Equations of motion were integrated via the velocity Verlet algorithm [40] and the time-step of the simulation was set at 0.5 fs. Periodic boundary conditions were applied inplane, in both the x and y directions to avoid the influence of the fixed walls, and the out-of-plane dimension of graphene was kept sufficiently large. Nose–Hoover thermostat [46,47] was applied to maintain the system at a constant temperature. Before calculating the heat flux and its autocorrelation function, we first equilibrated each structure for 0.2 ns (400,000 time steps) within the constant-(NVT) ensemble, with the temperature kept at T = 300 K. Following equilibration, we turned off the thermostat and the simulation was carried out in the constant-(NVE) ensemble for another 1million time steps (0.5 ns). Ultimately, our simulations were continued in the NVE ensemble to calculate the

At first, we calculated the thermal conductivity (TC) of an approximately square graphene sheet without nanomesh. TCs were computed by averaging the integral of heat flux autocorrelation function (HFACF) curves from twenty uncorrelated NVE ensembles with different initial conditions. Our result for the thermal conductivity, i.e., 2710 ± 140 W/mK, is consistent with that

Fig. 1. Schematic representation of the atomic structure of a graphene nanomesh and its supercell.

249

M. Yarifard et al. / Computational Materials Science 111 (2016) 247–251

Fig. 2. Graphene nanomeshes with different nanohole geometries but having the same area (same porosity).

Table 1 Comparison of thermal conductivities for different sizes of simulation box at room temperature. The porosity is 30% and the length of supercell is 5 nm. Number of supercell

System size (atoms)

Thermal conductivity (W/mK)

11 22 33 44

709 2733 5851 9939

31.1 ± 1.7 32.4 ± 1.3 30.7 ± 1.2 32.8 ± 1.7

Thermal Conductivity (W/mK)

400

Neck = 13 nm Neck = 7 nm Neck = 5 nm

350 300

(a)

250 200 150 100 50

5

10

15

20

25

30

35

40

45

Periodicity (nm) 400

Thermal Conductivity (W/mK)

obtained previously via the BTE method [50], and with experimental result [3] and with the EMD result [51]. Under the same simulation conditions, the thermal conductivity of the graphene nanomesh was calculated. Finite size effect appears when calculating thermal conductivity wherein the simulation box is not large enough. We calculated the thermal conductivity of GNMs with the porosity of 30%, which have different sizes, employing the EMD method at room temperature. As listed in Table 1, the values of the thermal conductivity are almost the same, within a 5% tolerance. Therefore, our simulation box was sufficiently large to overcome the finite size effect. In all our simulations, we used 4 supercells as the simulation box, i.e., 2 supercells in each direction. The periodicity and the neck width are indeed two important parameters for the design of a graphene nanomesh [28]. However, recent experimental [52] and theoretical [53] works on silicon nanomesh structures have found that periodicity and porosity are two dominant factors in determining the thermal conductivity. In this study, the periodicity (P), the porosity (u) and the neck width (N) are used as control parameters to clearly identify the key factors in determining the thermal conductivity. Fig. 3(a) represents the decreasing behavior and the final convergence of the thermal conductivity as a function of periodicity. Also, as can be seen, for a given periodicity, the thermal conductivity of GNM with a small neck width is low. To gain a better insight, we have shown the variations of the thermal conductivity with N/P (the ratio of neck width and periodicity). As is seen from Fig. 3(b), for N/P smaller than 0.4, the thermal conductivity is not affected by the periodicity of GNMs through the variation of N/P. Reduction in the thermal conductivity can occur

Neck = 13 nm Neck = 7 nm Neck = 5 nm

350 300 250

(b)

200 150 100 50

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Neck / Periodicity Fig. 3. Variations of thermal conductivity of a GNM with circular nanoholes (a) with periodicity (b) with neck/periodicity, for given neck widths.

M. Yarifard et al. / Computational Materials Science 111 (2016) 247–251

ðjporous =jbulk ÞRussell ¼ ð1  /2=3 Þ=ð1  /2=3 þ /Þ [53,56], indicating that the thermal conductivity should have the same value when the porosity is 43%. These models are only valid for structures in which the material dimensions are much larger than the phonon mean free path, and hence are not applicable to GNM, due to the strong size effect in nanostructures. Another way to modulate the phonon transport in the GNM is to vary the porosity. The dependence of thermal conductivity on the porosity is shown in Fig. 5. We selected the periodicity for all simulated GNMs to be 10 nm and varied the radii of nanoholes. Increasing the nanohole radius reduced the channels for phonon transport (neck width). Therefore, the increase of porosity resulted in the decrease of thermal conductivity. This rapid decay was due

400

Thermal Conductivity (W/mK)

because of two basic reasons. It may be due to a decrease in the phonon group velocity as a consequence of altering the phonon bonds (leading to phonon confinement) [54], or it can occur through modification of the scattering relaxation rates due to Umklapp process, and impurity and boundary scattering (although the phonon bond structure remains bulk-like). In this case, the radius of nanoholes increases as the periodicity enhances for a given neck width. Therefore, the decrease in the thermal conductivity is due to an increase in boundary scattering. The dependence of the thermal conductivity on the periodicity becomes negligible when N/P < 0.4, which is due to an enhancement of the ratio of the caging effect to the boundary scattering. It seems that the neck width and the size of nanoholes represent the caging effect and the boundary scattering, respectively. Tang et al. [52] have reported on an experimental investigation of thermal conductivity of holed thin silicon films. They found that by reducing the periodicity for a given porosity, the thermal conductivity was reduced significantly and the thermoelectric performance of holed silicon was better than the silicon nanowire systems so that nanoporous semiconductors may be highly attractive materials in thermoelectric applications [55]. Besides the periodicity effect for a given neck width, we also investigated the periodicity effect for a given porosity. Fig. 4 displays the dependence of the thermal conductivity on the periodicity of the GNM with 43% porosity. We see that the thermal conductivity is an increasing function of the periodicity. These results are consistent with previous studies [52]. This dependency deviates from the classical Eucken model ðjporous =jbulk ÞEucken ¼ ð1  /Þ=ð1  /=2Þ and the Russell model

Periodicity =10 nm

300

200

100

0 0

15

30

45

60

Porosity (%) Fig. 5. Thermal conductivity of a GNM with circular nanoholes as a function of porosity.

240

Thermal Conductivity (W/mK)

250

Porosity 22% Porosity 34%

200

Periodicity= 15nm

4

160

4

6

(a)

8

120

O8

O

6

80

3

3 40 3

4

5

6

7

8

Neck Width (nm) 240

240

Thermal Conductivity (W/mK)

Thermal Conductivity (W/mK)

Porosity 22% Porosity 34%

Porosity = 43 %

210 180 150 120 90

Eucken model = 1270 W/mK

60

200

Periodicity= 15nm

4

(b) 120

6

O

8

O

80

8

6

3

3

Russell model = 1355 W/mK

30

4

160

40 6

12

18

24

30

36

42

48

54

60

66

72

Periodicity (nm) Fig. 4. Thermal conductivity of a GNM with circular nanoholes as a function of periodicity.

3

4

5

6

7

8

Neck Width (nm) Fig. 6. Dependence of thermal conductivity on the neck width for different nanohole geometries (a) in x-direction (b) in y-direction.

M. Yarifard et al. / Computational Materials Science 111 (2016) 247–251

to an enhancement of boundary scattering and the decrease in the neck width that caused the phonons to be trapped. Snyder et al. [16] demonstrated that the reduction in the thermal conductivity of a silicon thin film through the formation of a nanomesh structure was due to alteration in the phonon band structure. Also, a previous Monte Carlo simulation of phonon transport in 2D porous silicon with square pores suggested that when the spacing between the adjacent pores is smaller than the phonon mean free paths (MFP), phonons can be trapped behind the pore [57]. We propose that both the boundary scattering and the phonon confinement play a role in reducing the thermal conductivity of GNMs, but the influence of the neck width on thermal conductivity is more pronounced. In other words, it seems that if the neck width is small, the phonons become trapped in cages between the nanoholes and the group velocity of heat-carrying agents decreases. Another useful issue to investigate was the geometry of nanoholes. We computed the thermal conductivity of GNMs with different geometries of the nanohole but with the same area and periodicity. Two collections with different nanohole areas were investigated. As illustrated in Fig. 2, the nanomeshes with different nanohole geometries have a different neck width. Fig. 6 displays the variations of thermal conductivity with the neck width of GNMs for different nanohole geometries, but with the same porosity for these two samples. We labeled the holes with triangular, square, hexagonal, octagonal and circular geometries with numbers 3, 4, 6, 8 and the letter O respectively. As is seen from Fig. 6, the thermal conductivity is affected by the nanohole geometries through the variation of the neck width in both directions. Since the neck width of the triangular case has the smallest value, the thermal conductivity for this case was obtained to be lower than those of the others. Also, the thermal conductivity was an increasing function with the neck width, and eventually the square case reveals a highest thermal conductivity in both directions. The simulation results show that the role of the neck width on TC is more pronounced, and that it seems that the caging of the phonons could be the responsible mechanism. 4. Conclusions EMD-based simulation was employed to investigate the thermal conductivity of GNMs. The computed results demonstrate that the thermal conductivity is reduced by more than one order of magnitude compared to the pristine graphene. If the high electronic conduction of graphene can be preserved due to the periodic arrangement of holes as opposed to random porous systems, the reduction in thermal conduction could be utilized in thermoelectric applications since it leads to an enhanced thermoelectric figure of merit (ZT) [15]. Neck width effect on TC is more pronounced and it seems that a caging effect of phonons dominates the boundary scattering for N/P < 0.4. Also, The influence of nanoholes geometry on GNM was explored and we found that a GNM with triangular nanoholes has the smallest thermal conductivity for a given porosity. Therefore, porosity, periodicity, neck width and the geometry of nanoholes can be used to tune the thermal conductivity of GNM for thermal managements. Acknowledgment M.Y. thanks A. Rajabpour from Imam Khomeini International University at Qazvin, and H. Alizadeh and F. Yousefi for useful discussions.

251

References [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666–669. [2] C. Soldano, A. Mahmood, E. Dujardin, Carbon 48 (2010) 2127–2150. [3] S. Ghosh, I. Calizo, D. Teweldebrhan, E.P. Pokatilov, D.L. Nika, A.A. Balandin, W. Bao, F. Miao, C.N. Lau, Appl. Phys. Lett. 92 (2008) 151911. [4] M.H. Shih, L.J. Li, Y.C. Yang, H.Y. Chou, C.T. Lin, C.Y. Su, ACS Nano 7 (2013) 10818–10824. [5] Y.M. Lin, C. Dimitrakopoulos, K.A. Jenkins, D.B. Farmer, H.Y. Chiu, A. Grill, P. Avouris, Science 327 (2010) 662. [6] Q. He, H.G. Sudibya, Z. Yin, S. Wu, H. Li, F. Boey, W. Huang, P. Chen, H. Zhang, ACS Nano 4 (2010) 3201–3208. [7] T.H. Han, Y.K. Huang, A.T. Tan, V.P. Dravid, J. Huang, J. Am. Chem. Soc. 133 (2011) 15264–15267. [8] Z. Yin, S. Wu, X. Zhou, X. Huang, Q. Zhang, F. Boey, H. Zhang, Small 6 (2010) 307–312. [9] S. Liu, T.H. Zeng, M. Hofmann, E. Burcombe, J. Wei, R. Jiang, J. Kong, Y. Chen, ACS Nano 5 (2011) 6971–6980. [10] X. Zhao, C.M. Hayner, M.C. Kung, H.H. Kung, ACS Nano 5 (2011) 8739–8749. [11] T.M. Tritt, Thermal conductivity: theory, properties and applications. New York, 2004. [12] C.W. Chang, D. Okawa, A. Majumdar, A. Zettl, Science 314 (2006) 1121–1124. [13] J. Hu, X. Ruan, Y.P. Chen, Nano Lett. 9 (2009) 2730–2735. [14] A. Majumdar, Science 303 (2004) 777–778. [15] H. Sevinçli, G. Cuniberti, Phys. Rev. B 81 (2010) 113401. [16] J. Yu, S. Mitrovic, D. Tham, J. Varghese, J.R. Heath, Nat. Nanotechnol. 5 (2010) 718–721. [17] A.A. Balandin, Nat. Mater. 10 (2011) 569–581. [18] H. Zhang, G. Lee, K. Cho, Phys. Rev. B 84 (2011) 115460. [19] T.Y. Ng, J.J. Yeo, Z.S. Liu, Carbon 50 (2012) 4887–4893. [20] D. Yanga, F. Maa, Y. Suna, T. Hua, K. Xu, Appl. Surf. Sci. 258 (2012) 9926–9931. [21] S.K. Chien, Y.T. Yang, C.K. Chen, Carbon 50 (2012) 421–428. [22] N. Wei, L. Xu, H.Q. Wang, J.C. Zheng, Nanotechnology 22 (2011) 105705. [23] J. Hu, S. Schiffli, A. Vallabhaneni, X. Ruan, Y.P. Chen, Appl. Phys. Lett. 97 (2010) 133107. [24] B. Mortazavi, A. Rajabpour, S. Ahzi, Y. Remond, S.M.V. Allaei, Solid State Commun. 152 (2012) 261–264. [25] Y. Wang, B. Qiu, X. Ruan, Appl. Phys. Lett. 101 (2012) 013101. [26] A. Cao, J. Qu, J. Appl. Phys. 111 (2012) 053529. [27] H.Y. Cao, H. Xiang, X.G. Gong, Solid State Commun. 152 (2012) 1807–1810. [28] J. Bai, X. Zhong, S. Jiang, Y. Huang, X. Duan, Nat. Nanotechnol. 5 (2010) 190– 194. [29] L. Jiang, Z. Fan, Nanoscale 6 (2014) 1922–1945. [30] V.H. Nguyen, F. Mazzamuto, J.S. Martin, A. Bournel, P. Dollfus, Nanotechnology 23 (2012) 065201. [31] A. Esfandiar, N.J. Keybert, E.N. Dattoli, G.H. Han, M.B. Lerner, O. Akhavan, A. Irajizad, A.T.C. Johnson, Appl. Phys. Lett. 103 (2013) 183110. [32] O. Akhavan, E. Ghaderi, Small 9 (2013) 3593–3601. [33] S. Blankenburg, M. Bieri, R. Fasel, K. Mullen, C.A. Pignedoli, D. Passerone, Small 6 (2010) 2266–2271. [34] S. Jungthawan, P. Reunchan, S. Limpijumnong, Carbon 54 (2013) 359–364. [35] A. Du, Z. Zhu, S.C. Smith, J. Am. Chem. Soc. 132 (2010) 2876–2877. [36] L. Hu, D. Maroudas, J. Appl. Phys. 116 (2014) 184304. [37] C. Carpenter, A.M. Christmann, L. Hu, I. Fampiou, A.R. Muniz, A. Ramasubramaniam, D. Maroudas, Appl. Phys. Lett. 104 (2014) 141911. [38] N.C.B. Mosterio, A.F. Fonseca, Phys. Rev. B. 89 (2014) 195437. [39] H. Sahin, S. Ciraci, Phys. Rev. B. 84 (2011) 035452. [40] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [41] J.M. Haile, Molecular Dynamics Simulation, Wiley-Interscience Publication, New York, 1992. [42] L. Lindsay, D.A. Broido, Phys. Rev. B 81 (2010) 205441. [43] P.K. Schelling, S.R. Phillpot, P. Keblinski, Phys. Rev. B 65 (2002) 144306. [44] S. Plimpton, J. Comp. Phys. 117 (1995) 1–19. [45] See http://lammps.sandia.gov for more information about LAMMPS. [46] S. Nose, J. Chem. Phys. 81 (1984) 511–519. [47] W.G. Hoover, Phys. Rev. A 31 (1985) 1695–1697. [48] D.W. Bullett, J. Phys. C: Solid State Phys. 8 (1975) 2707–2714. [49] J. Chen, G. Zhang, B. Li, Phys. Lett. A 374 (2010) 2392–2396. [50] D.L. Nika, E.P. Pokatilov, A.S. Askerov, A.A. Balandin, Phys. Rev. B 79 (2009) 155413. [51] H.K. Liu, Y. Lin, S.N. Luo, J. Phys. Chem. C 118 (2014) 24797–24802. [52] J. Tang, H.T. Wang, D.H. Lee, M. Fardy, Z. Huo, T.P. Russell, P. Yang, Nano Lett. 10 (2010) 4279–4283. [53] L. Yang, N. Yang, B. Li, Nano Lett. 14 (2014) 1734–1738. [54] A. Balandin, K.L. Wang, Phys. Rev. B 58 (1998) 1544–1549. [55] J.H. Lee, G.A. Galli, J.C. Grossman, Nano Lett. 8 (2008) 3750–3754. [56] H.W. Russell, J. Am. Ceram. Soc. 18 (1935) 1–5. [57] Q. Hao, G. Chen, M.S. Jeng, J. Appl. Phys. 106 (2009) 114321.