Abstract
Thermal transport properties of graphene with nanosized constrictions are investigated using nonequilibrium molecular dynamics simulations. The results show that the nanosized constrictions have a significant influence on the thermal transport properties of graphene. The thermal resistance of the nanosized constrictions is on the order of 10^{7} to 10^{9} K/W at 150 K, which reduces the thermal conductivity by 7.7% to 90.4%. It is also found that the constriction resistance is inversely proportional to the width of the constriction and independent of the heat current. Moreover, we developed an analytical model for the ballistic thermal resistance of the nanosized constrictions in twodimensional nanosystems. The theoretical prediction agrees well with the simulation results in this paper, which suggests that the thermal transport across the nanosized constrictions in twodimensional nanosystems is ballistic in nature.
PACS
65.80.CK; 61.48.Gh; 63.20.kp; 31.15.xv
Keywords:
Graphene; Ballistic resistance; Nanosized constriction; Molecular dynamics simulationBackground
Graphene is a twodimensional (2D) material formed of the honeycomb lattice of sp^{2}bonded carbon atoms. The strong bonding and perfect lattice structure give its unique thermal properties [13]. As Balandin et al. [1,2] demonstrated, the thermal conductivity of graphene is up to 5,400 W/(m · K), which makes it one of the most promising base materials for nextgeneration electronics and thermal management [26]. Additionally, compared with other highconductivity materials, such as carbon nanotubes [79], graphene is much easier to be fashioned into a broad range of shapes. Such flexibility makes possible the utilization of graphene.
Usually, limited by the synthesis and fabrication procedure, graphene inevitably has a variety of defects, such as vacancies, StoneWales defects, and impurities [10,11]. Many scholars have demonstrated that these defects are obstacles to heat transfer and create additional sources of phonon scattering in graphene [1216], especially when the characteristic dimension is less than the phonon mean free path (approximately 775 nm at room temperature) [2]. Hao et al. [13] performed molecular dynamics (MD) simulations on defected graphene sheets. They observed that the increasing defect concentration dramatically reduces the thermal conductivity of graphene. Chien et al. [14] investigated the effect of impurity atoms in graphene and found a rapid drop in thermal conductivity, where hydrogen coverage down to as little as 2.5% of the carbon atoms reduces the thermal conductivity by about 40%. So we can conclude that the thermal transport properties of graphene are very sensitive to its own structures. Besides these defects, the structural configuration is another important but less studied factor impacting the thermal properties, and thus, it can affect the lifetime and reliability of the graphenebased nanodevices further because these devices have more complex shapes in engineering situations. Therefore, from a practical point of view, the investigation on how to predict or tune the thermal transport properties of graphene with a variety of shapes is especially useful for thermal management.
Recently, Xu et al. [17] investigated the transport properties of various graphene junctions and quantum dots using nonequilibrium Green's function method and found that the thermal conductance is insensitive to the detailed structure of the contact region but substantially limited by the narrowest part of the system. Huang et al. [18] constructed a sandwich structure with atomistic Green's function method, where two semiinfinite graphene sheets are bridged by a graphene nanoribbon (GNR). They mainly focused on the phonon transport behavior in GNR and observed that the thermal conductance increases with the width of GNR at fixed length and decreases with GNR length at fixed width.
This paper presents the effect of the nanosized constrictions on the thermal transport properties of graphene studied by the nonequilibrium molecular dynamics (NEMD) simulations. We calculate the thermal transport properties of graphene with those constrictions, and the effects of the heat current and the width of the constriction were explored in detail. Further, based on the phonon dynamics theory, we develop an analytical model for the ballistic resistance of the nanosized constrictions in twodimensional nanosystems, which agrees well with the simulation results in this paper.
Methods
Here, we employed the NEMD method [1924] to simulate the thermal transport in graphene. The simulated system with constrictions is illustrated in Figure 1, which is originally an 18.2nmlong and 11.9nmwide rectangular graphene sheet with zigzag long edges. Fixed boundary conditions are used at the outmost layers of each end along the length direction, i.e., the green atoms in Figure 1, to prevent spurious global rotation and translation of the graphene. Free boundary conditions are used along the width direction. As depicted in Figure 1, in the middle of the system, three nanosized constrictions are constructed by introducing four linear vacancy defects into the graphene sheet, so that the thermal transport is possible only through the small area in contact. These constrictions are in the same size and distribute uniformly along the width direction. As shown in Figure 1b, the width of one constriction is w = (w_{1} + w_{2})/2 and the total cross section area of three constrictions is A = 3wδ, in which δ = 0.335 nm is the thickness of the graphene sheet [3,25].
Figure 1. Schematic of molecular dynamics simulation. (a) Simulation system including a hightemperature slab (red) and a lowtemperature slab (blue) with fixed boundaries (green). (b) Detailed structure of the constriction.
In the MD simulations, the bondorder potential presented by Brenner [26] is used to describe the carboncarbon bonding interactions,
where E_{b} is the total potential energy, V_{R} and V_{A} are the pairadditive repulsive and attractive potential terms, respectively, f(r_{ij}) is the truncation function that explicitly restricts the potential to nearest neighbors, and b_{ij} is the manybody interaction parameter. The atomic motion is integrated by a leapfrog scheme with a fixed time step of 0.5 fs. Each simulation case runs for 1 ns to reach a steady state, and then for 1.5 ns to average the temperature profile and heat current over time. During the simulation, the mean temperature of all cases is set at 150 K, which is maintained by the NoséHoover thermostat method [27]. The heat current is generated by exchanging the velocity vector of one atom in the hightemperature slab (the red part) and another in the lowtemperature slab (the blue part) constantly. This method was developed by MüllerPlathe [28], and it can keep the total energy and momentum of the system conserved. The heat current is defined as
in which m is the atomic mass of carbon, v_{h} is the velocity of the hottest atom in the lowtemperature slab, v_{c} is the velocity of the coldest atom in the hightemperature slab, and t is the statistical time. Specifically, by comparing the actual heat current with the preset heat current, we can adjust the frequency of the velocity exchange in real time and achieve that preset heat current finally. After reaching steady state, the system is equally divided into 50 slabs along the length direction. And the local instantaneous temperature for each slab is defined through the averaged kinetic energy according to the energy equipartition theorem as
where N is the number of atoms per slab, k_{B} is the Boltzmann constant, and P_{i} is the momentum of the ith atom.
Results and discussion
Nine graphene sheets with differentsized constrictions are simulated in this paper, and the corresponding pristine one is also designed for comparison. The constriction widths of nine cases are 0.216, 0.648, 1.08, 1.512, 1.944, 2.376, 2.808, 3.24, and 3.672 nm, respectively. And four heat currents (i.e., J = 0.2097, 0.3146, 0.4195, and 0.5243 μW) are performed for all the cases.
The typical temperature profile of the graphene with nanosized constrictions is shown in Figure 2. As mentioned before, we produce an energy transfer from the sink region to the source region by exchanging the velocities. Therefore, several additional phonon modes are excited, which leads to the temperature jumps near the high and lowtemperature slabs [29]. Between those slabs and constrictions, the temperature distribution is linear, but not completely symmetrical. Specifically, on the left side of the system, the mean temperature is 175 K and the thermal conductivity calculated by the Fourier law is 110 W/(m · K), while on the right side, the mean temperature is 125 K and a higher thermal conductivity, 133 W/(m · K), is obtained. The asymmetry shows the obvious temperature dependence of the thermal conductivity of graphene, which is consistent with the results confirmed by Balandin et al. on the aspects of firstprinciple calculations and experiments [1,12]. Besides, in the following, we will mainly focus on the big temperature jump ∆T at the constriction as shown in Figure 2, which indicates that energy is blocked when passing through the constriction and thus an additional thermal resistance is introduced.
Figure 2. Typical temperature profile. The temperature profile is obtained by injecting the heat current of 0.5243 μW. The inset shows the corresponding simulation system with the constriction width of 1.512 nm.
The temperature profiles of the systems with differentsized constrictions, under different heat current, are shown in Figure 3. And the insets show the dependence of the temperature jump ∆T extracted from those temperature profiles on the heat current. As shown in Figure 3, with the heat current increasing, the temperature jump approximately increases linearly, which indicates that the thermal resistance at the constrictions is an intrinsic property of the system and it is independent of the heat current, while for different systems, with a fixed heat current, the temperature jump varies with the constriction width. When the width is 1.08 nm, the temperature jump spans the range 25.5 to 63 K. But when the width is 1.512 nm, the range is from 18 to 42 K, onethird lower than the former. This thermal transport behavior is distinctly different from that of the bulk material, which is independent of the size, and indicates that the thermal resistance of constriction in graphene has obvious size effects.
Figure 3. Temperature profiles versus heat current. (a, b) From different systems with the constriction widths of 1.08 and 1.512 nm, respectively. The insets show the temperature jump extracted from the temperature profiles versus the heat current.
Similar to the interfacial thermal resistance, i.e., Kapitza resistance, the thermal resistance R at the constrictions can be defined as
where J and ∆T, respectively, correspond to the heat current across the constrictions and the associated temperature jump (as shown in Figure 2). In order to reduce the error, in this paper, the constriction resistance R is calculated by fitting the curve between the temperature jump and the heat current. The results are shown in Figure 4, where w is the width of one constriction, with larger w meaning weaker strength of the constriction. The results show that the nanosized constriction resistance is on the order of 10^{7} to 10^{9} K/W. And as mentioned before, the constriction resistance has an obvious size effect, which decreases from 4.505 × 10^{8} to 9.897 × 10^{6} K/W with the increasing width, and it is almost inversely proportional to the width of the constrictions.
Figure 4. Constriction resistance versus width of constriction. The dots are MD results and the curve is the theoretical prediction given by Equation 9.
To quantitatively describe the effect of the nanosized constrictions on thermal transport properties, we introduce a dimensionless parameter: the thermal conductance ratio η = σ/σ_{0}, where σ and σ_{0} are the thermal conductance of the graphene with constrictions and that of the corresponding pristine graphene, respectively. Figure 5 shows the dependence of the thermal conductance ratio on the width. As shown, varioussized constrictions have a significant influence on the thermal conductance of graphene and the thermal conductance is reduced by 7.7% to 90.4%. Thus, we can conclude that it is quite feasible to tune the thermal conductance of graphene over a wide range by introducing the nanosized constriction or controlling the configuration of the embedded extended defect in graphene.
Figure 5. Thermal conductance ratio versus width of constriction. The inset is the corresponding pristine graphene.
Recently, some modelbased analyses on the constriction resistance have been carried out [3033]. The models mainly involve the following three parameters: the phonon mean free path (l), the characteristic size of the constriction (a), and the dominant phonon wavelength (λ_{d}). In the completely diffusive regime when a is much larger than l, the diffusive constriction resistance (R_{d}) is given by the Maxwell constriction resistance model [30]:
where κ denotes the thermal conductivity. But in the other limit, that is, a < < l, phonon transport across the constriction is ballistic. The heat current in the ballistic regime [31,32] can be written as
where ω is the frequency of phonons, ω_{m} is the maximum frequency, ℏ is the reduced Planck constant, is the occupation of phonons given by the BoseEinstein distribution, D(ω) is the phonon density of states, v_{g}(ω) is the phonon group velocity, τ(ω,θ) is the transmissivity of phonons, θ is the polar angle, and φ is the azimuthal angle. What is more, in the ballistic limit, two limiting cases of phonon transmission behavior are further discussed, which is differentiated depending on the characteristic size of the constriction (a) relative to the dominant phonon wavelength λ_{d}. If a is much larger than λ_{d}, which is the geometric scattering limit, the transmissivity of phonons is described as τ(ω,θ) = cosθ. If a is near or smaller than λ_{d}, which is the Rayleigh scattering limit, the effect of the wave diffraction should be considered and the calculation of the transmissivity is more complex [33]. It can be seen that the theoretical modeling of the constriction resistance is based on the threedimensional (3D) system so far. But for graphene, a 2D material, it is invalid.
In this paper, the width of one constriction in graphene is 0.216 ~ 3.672 nm, which is much smaller than the phonon mean free path of graphene (approximately 775 nm) with 2 orders of magnitude. Therefore, the thermal transport at the constrictions is in the ballistic regime. In analogy to the 3D ballistic model, the heat current for 2D nanosystems can be described as
where the dominant phonon wavelength is λ_{d} ≈ 2.3hv_{g}/(k_{B}T) [33], in which h is the Planck constant. We assume that the phonon group velocity (v_{g}) is independent of phonon modes and frequency. Then we get λ_{d} = 12.84 nm by substituting the phonon group velocity v_{g} = 17.45 km/s (the average of v_{LA} = 21.3 km/s for the LA mode and v_{TA} = 13.6 km/s for the TA mode in graphene [12]). Therefore, the transmissivity of phonons is τ(ω,θ) = cosθ, and Equation 7 can be simplified to
where U is the internal energy per unit volume. Thus, the ballistic constriction resistance of the 2D nanosystems is
From Equation 9, the ballistic constriction resistance is inversely proportional to the cross section area (A), i.e., the width of the constriction (w), which is consistent with the conclusion of MD. And the predicted results, obtained by substituting c_{v} = 6.81 × 10^{5} J/(m^{3} · K) [34] and v_{g} = 17.45 km/s into Equation 9, are compared quantitatively with MD results in Figure 4. It can be seen that the present model predicts well the thermal resistance of the constriction in graphene, which suggests that thermal transport across the nanosized constrictions in 2D nanosystems is ballistic in nature.
Conclusions
Graphene has shown great potential for the applications in highefficiency thermal management and nanoelectronics due to its exceptional thermal properties in the past few years. Understanding the underlying mechanism of controlling the thermal properties of various structures is of considerable interest. In this paper, systems of rectangular graphene sheets with various nanosized constrictions are constructed by embedding linear vacancy defects and the thermal transport properties are investigated by using nonequilibrium molecular dynamics method. The results show that the nanosized constriction has a significant influence on the thermal transport properties of graphene. And the constriction resistance is on the order of 10^{7} to 10^{9} K/W at 150 K, which reduces the thermal conductivity by 7.7% to 90.4%. Besides, the constriction resistance is inversely proportional to the constriction width and independent of the heat current. These findings indicate that the desired thermal conduction can be achieved via nanosized constrictions. Moreover, we develop a ballistic constriction resistance model for 2D nanosystems, which corresponds to the case when the mean free path of phonon is much larger than the characteristic dimension of the constriction. The predicted values of this model agree well with the simulation results in this paper, which suggests that the thermal transport across nanosized constrictions in 2D nanosystems is ballistic in nature.
Abbreviations
2D: twodimensional; 3D: threedimensional; GNR: graphene nanoribbon; MD: molecular dynamics; NEMD: nonequilibrium molecular dynamics.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
BYC conceived of the study; participated in its design, coordination, and analyses; and revised the manuscript critically for important intellectual content. WJY carried out the molecular dynamics simulations, interpreted the results, and drafted the manuscript. HMY and BMC performed the data analyses and edited the manuscript critically. All authors discussed the results and read and approved the final manuscript.
Acknowledgements
This work was supported by the National Natural Science Foundation of China (Grant Nos. 51322603, 51136001, and 51356001), Science Fund for Creative Research Groups (No. 51321002), the Program for New Century Excellent Talents in University, Tsinghua University Initiative Scientific Research Program, the Tsinghua National Laboratory for Information Science and Technology of China, and the Foundation of Key Laboratory of Renewable Energy Utilization Technologies in Buildings of the National Education Ministry in Shandong Jianzhu University (No. KF201301).
References

Balandin AA, Ghosh S, Bao W, Calizo I, Teweldebrhan D, Miao F, Lau CN: Superior thermal conductivity of singlelayer graphene.
Nano Lett 2008, 8:902907. PubMed Abstract  Publisher Full Text

Ghosh S, Calizo I, Teweldebrhan D, Pokatilov EP, Nika DL, Balandin AA, Bao W, Miao F, Lau CN: Extremely high thermal conductivity of graphene: prospects for thermal management applications in nanoelectronic circuits.
Appl Phys Lett 2008, 92:15191113. Publisher Full Text

Pop E, Varshney V, Roy AK: Thermal properties of graphene: fundamentals and applications.
MRS Bull 2012, 37:12731281. Publisher Full Text

Novoselov KS, Geim AK, Morozov SV, Jiang D, Zhang Y, Dubonos SV, Grigorieva IV, Firsov AA: Electric field effect in atomically thin carbon films.
Science 2004, 306:666669. PubMed Abstract  Publisher Full Text

Geim AK, Kim P: Carbon wonderland.
Sci Am 2008, 298:9097. PubMed Abstract

Soldano C, Mahmood A, Dujardin E: Production, properties and potential of graphene.
Carbon 2010, 48:21272150. Publisher Full Text

Fujii M, Zhang X, Xie H, Ago H, Takahashi K, Ikuta T, Abe H, Shimizu T: Measuring the thermal conductivity of a single carbon nanotube.
Phys Rev Lett 2005, 95:06550214. PubMed Abstract  Publisher Full Text

Pop E, Mann D, Wang Q, Goodson K, Dai H: Thermal conductance of an individual singlewall carbon nanotube above room temperature.
Nano Lett 2006, 6:96100. PubMed Abstract  Publisher Full Text

Choi TY, Poulikakos D, Tharian J, Sennhauser U: Measurement of the thermal conductivity of individual carbon nanotubes by the fourpoint threeω method.
Nano Lett 2006, 6:15891593. PubMed Abstract  Publisher Full Text

Hashimoto A, Suenaga K, Gloter A, Urita K, Iijima S: Direct evidence for atomic defects in graphene layers.
Nature 2004, 430:870873. PubMed Abstract  Publisher Full Text

Lee GD, Wang CZ, Yoon E, Hwang NM, Kim DY, Ho KM: Diffusion, coalescence, and reconstruction of vacancy defects in graphene layers.
Phys Rev Lett 2005, 95:20550114. PubMed Abstract  Publisher Full Text

Nika DL, Pokatilov EP, Askerov AS, Balandin AA: Phonon thermal conduction in graphene: role of Umklapp and edge roughness scattering.

Hao F, Fang D, Xu Z: Mechanical and thermal transport properties of graphene with defects.
Appl Phys Lett 2011, 99:04190113. Publisher Full Text

Chien S, Yang Y, Chen C: Influence of hydrogen functionalization on thermal conductivity of graphene: nonequilibrium molecular dynamics simulations.
Appl Phys Lett 2011, 98:03310713. Publisher Full Text

Yang P, Wang XL, Li P, Wang H, Zhang LQ, Xie FW: The effect of doped nitrogen and vacancy on thermal conductivity of graphenenanoribbon from nonequilibrium molecular dynamics.
Acta Phys Sin 2012, 61:07650118.
in Chinese

Yao HF, Xie YE, Tao O, Chen YP: Thermal transport of graphene nanoribbons embedding linear defects.
Acta Phys Sin 2013, 62:06810217.
in Chinese

Xu Y, Chen X, Wang JS, Gu BL, Duan W: Thermal transport in graphene junctions and quantum dots.

Huang Z, Fisher TS, Murthy JY: Simulation of thermal conductance across dimensionally mismatched graphene interfaces.
J Appl Phys 2010, 108:11431017. Publisher Full Text

Ye ZQ, Cao BY, Guo ZY: High and anisotropic thermal conductivity of bodycentered tetragonal C_{4} calculated using molecular dynamics.

Hu GJ, Cao BY: Thermal resistance between crossed carbon nanotubes: molecular dynamics simulations and analytical modeling.
J Appl Phys 2013, 114:22430818. Publisher Full Text

Li YW, Cao BY: A uniform sourceandsink scheme for calculating thermal conductivity by nonequilibrium molecular dynamics.
J Chem Phys 2010, 133:02410615. PubMed Abstract  Publisher Full Text

Hu GJ, Cao BY: Molecular dynamics simulations of heat conduction in multiwalled carbon nanotubes.
Mol Simulat 2012, 38:823829. Publisher Full Text

Cao BY, Kong J, Xu Y, Yung K, Cai A: Polymer nanowire arrays with high thermal conductivity and superhydrophobicity fabricated by a nanomolding technique.
Heat Transfer Eng 2013, 34:131139. Publisher Full Text

Yao WJ, Cao BY: Thermal wave propagation in graphene studied by molecular dynamics simulations.

Hu J, Ruan X, Chen YP: Thermal conductivity and thermal rectification in graphene nanoribbons: a molecular dynamics study.
Nano Lett 2009, 9:27302735. PubMed Abstract  Publisher Full Text

Brenner DW: Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films.
Phys Rev B 1990, 42:94589471. Publisher Full Text

Hoover WG: Canonical dynamics: equilibrium phasespace distributions.
Phys Rev A 1985, 31:16951697. PubMed Abstract  Publisher Full Text

MüllerPlathe F: A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity.
J Chem Phys 1997, 106:60826085. Publisher Full Text

Jiang JW, Chen J, Wang JS, Li BW: Edge states induce boundary temperature jump in molecular dynamics simulation of heat conduction.

Cooper MG, Mikic BB, Yovanovish MM: Thermal contact conductance.
Int J Heat Mass Transfer 1969, 12:279300. Publisher Full Text

Prasher R: Predicting the thermal resistance of nanosized constrictions.
Nano Lett 2005, 5:21552159. PubMed Abstract  Publisher Full Text

Prasher R: Ultralow thermal conductivity of a packed bed of crystalline nanoparticles: a theoretical study.

Prasher R, Tong T, Majumdar A: Diffractionlimited phonon thermal conductance of nanoconstrictions.
Appl Phys Lett 2007, 91:14311913. Publisher Full Text

Mounet N, Marzari N: Firstprinciples determination of the structural, vibrational and thermodynamic properties of diamond, graphite, and derivatives.