Abstract

Single-layer graphene is so flexible that its flexural mode (also called the ZA mode, bending mode, or out-of-plane transverse acoustic mode) is important for its thermal and mechanical properties. Accordingly, this review focuses on exploring the relationship between the flexural mode and thermal and mechanical properties of graphene. We first survey the lattice dynamic properties of the flexural mode, where the rigid translational and rotational invariances play a crucial role. After that, we outline contributions from the flexural mode in four different physical properties or phenomena of graphene – its thermal conductivity, thermal expansion, Young’s modulus, and nanomechanical resonance. We explain how graphene’s superior thermal conductivity is mainly due to its three acoustic phonon modes at room temperature, including the flexural mode. Its coefficient of thermal expansion is negative in a wide temperature range resulting from the particular vibration morphology of the flexural mode. We then describe how the Young’s modulus of graphene can be extracted from its thermal fluctuations, which are dominated by the flexural mode. Finally, we discuss the effects of the flexural mode on graphene nanomechanical resonators, while also discussing how the essential properties of the resonators, including mass sensitivity and quality factor, can be enhanced.

This review focuses on the connection between the fundamental lattice dynamics and the thermal and mechanical properties of the unique two-dimensional material graphene. The lattice dynamical properties, which comprise the phonon spectrum or phonon modes, give fundamental information regarding the atomic interaction within the material. The framework for lattice dynamics was established by Born in the 1920s, and further developed by Debye, Einstein, Mott and others in the following decades. An interesting point is that while lattice dynamics is an atomic scale theory, direct connections can be made to macroscopic phenomena and properties. Readers are referred to the book by Born and Huang for a comprehensive description and discussion regarding lattice dynamic theory.Born and Huang (1954)

We focus on the flexural mode, which is a characteristic vibration mode in solid plates or rods.Landau and Lifshitz (1995) Elastic waves in solid plates are guided by the two outer surfaces, at which the component of stress is zero in the perpendicular direction. Kichhoff formulated the exact equation and boundary conditions for the flexural vibration of a plate in 1850.Kirchhoff (1850) A practical theoretical model of the flexural motion was first developed by Rayleigh in 1885.Rayleigh (1885) In 1917, Lamb completely solved the surface vibration problem.Lamb (1917) It was found that the free surface boundary condition restricts the elastic waves in the plate into two infinite sets of Lamb waves. The motion of one set of the Lamb wave is symmetrical about the midplane of the plate, while the other is anti-symmetric about the midplane. The zero-order (lowest-frequency) wave from the antisymmetric Lamb set is the flexural mode, whose defining characteristic is its parabolic dispersion curve. In this sense, flexural mode is one particular form of the Lamb wave in solid plates. Because the flexural mode has the lowest frequency among all Lamb waves, it is the easiest to be excited and carries most of the vibrational energy.

In microscopic lattice dynamics theory, the surface phonon mode is well known, and the long wave limits of the flexure phonon mode is the acoustic flexural wave in elastic mechanics. For three-dimensional bulk solid with volume L3, the ratio of surface phonon mode number to total phonon mode number is proportional to 1/L. Hence, in large piece of bulk materials, surface phonon modes (including the flexural mode) are not so important, except for some specific surface-related topics such as surface reconstruction, surface adsorption, and surface chemistry, etc. However, nanomaterials have large surface to volume ratio, so surface modes play an important role. As an extreme case, all atoms in graphene are exposed on the surface, so all phonon modes are essentially “surface modes”, and the flexural mode occupies one sixth of all phonon modes. That is the origin for the importance of the flexural mode in graphene.

The understanding and study of flexural modes in other novel materials have intensified significantly in recent years due to the recent isolation of the thinnest possible two-dimensional, one-atom-thick “plate”, graphene. Graphene is the best-known one-atom-thick two-dimensional material, which earned a Nobel prize in physics for Novoselov and Geim in 2010.Geim and Novoselov (2007) Accordingly, substantial effort has been expended in achieving a basic understanding of the lattice dynamical properties of single-layer graphene due to its status as the thinnest possible two-dimensional material. A key outcome of its low-dimensional structure is the existence of a flexural mode (also named ZA mode, bending mode, or out-of-plane transverse acoustic mode) in graphene. As the long wave flexural mode has the lowest frequency among all phonon modes, it is the easiest to be excited in graphene.

A complete understanding of the effect that the flexural mode has on the thermal and mechanical properties of graphene will be essential for many of the key applications graphene has been envisioned for. For example, graphene has been touted as the next-generation replacement for silicon in integrated circuits, though its bandgap is zero unless mechanical strain, doping, or finite width nanoribbons are created.Novoselov et al. (2005a) For these electronic applications, it is important to efficiently remove heat from the transistor during its high speed operation. Graphene possesses a superior thermal conductivity, which is crucial to prevent graphene based transistors from suffering thermally-induced malfunctioning during operation.Balandin et al. (2008); Nika et al. (2009a, b); Balandin (2011); Jiang et al. (2009a) It is now clear that the high thermal conductivity of graphene is contributed by its three acoustic modes at room temperature, including flexural mode.

The flexural mode also controls the behavior and properties of graphene nanoelectromechanical systems (NEMS) and nanoresonators, which have been proposed for various sensing applications, in particular ultrasensitive mass sensing and detection. However, the performance of these NEMS, and in particular its quality factor and energy dissipation, is strongly controlled by the flexural mode. Therefore, understanding how the limits imposed by the flexural mode as well as methods to circumvent these flexural mode-related loss mechanisms will be essential to enabling graphene NEMS devices.

From the above, it is clear that the flexural mode plays a critical role in governing the thermal and mechanical properties of graphene that will be critical for the success of graphene as a commercially-viable material. Due to the ongoing interest in graphene, it is an appropriate time to review the topic of the flexural mode in graphene, and its impact on graphene’s physical properties.

ii.1 Introduction

The phonon dispersion of graphene is calculated through the diagonalization of its dynamical matrix, which is a 6×6 complex matrix, due to two inequivalent carbon atoms in its primitive unit cell. The dynamical matrix is constructed based on the interatomic interaction and the space group symmetry of the system.Born and Huang (1954) The eigenvalue solution of the dynamical matrix gives the frequency and the eigen vector (vibration morphology) for all phonon modes in the system. The vibration morphology of the phonon mode actually represents a particular symmetric mechanical vibration of the graphene. Phonon modes can be used to decompose a general movement of the graphene. In this sense, phonon modes are usually called normal modes.

The lattice dynamics provides fundamental information for mechanical and thermal properties in graphene. For instance, the flexural mode is characteristic by its parabolic spectrum, which is related to its superior thermal conductivity.Balandin et al. (2008) The special vibration morphology of the flexural mode is the origin for the strong thermal contraction effect in graphene in a wide temperature range.Bao et al. (2009) The flexural mode’s long life time results in a high quality (Q) factor of the graphene nanomechanical resonator (GNMR).Jensen et al. (2008) These few examples demonstrate the importance of the flexural mode in graphene. Hence, our first task in present review is to recall the computation of the phonon dispersion for graphene.

ii.2 Lattice Dynamics for Flexural Mode

Bloch’s Theorem

In this subsection, we take the honeycomb lattice structure of graphene as an example to explain the application of the Bloch’s theorem. Graphene has a honeycomb lattice structure as shown in Fig. 1. Its symmetry is described by the D6h point symmetry group.Tang et al. (2011) According to its space group, the whole honeycomb lattice can be obtained by repeating a rhombus primitive unit cell, with two bases →a1 and →a2. The lattice constantSaito et al. (1998)|→ai|=2.46 Å. Each cell is denoted by a pair of integers (l1, l2). Carbon atoms are indexed by (l1, l2, s), where s=A or B are the two inequivalent carbon atoms in each unit cell. The locations of atoms A and B in the unit cell (0, 0) are →τA=(→a1+→a2)/3 and →τB=2(→a1+→a2)/3. The position of an arbitrary carbon atom is determined by →r(l1,l2,s)=→Rl1l2+→τs. The lattice vector is →Rl1l2=l1→a1+l2→a2.

The structure of an infinite graphene sheet is unchanged after it is shifted for a lattice vector →Rl1l2. That is graphene has lots of translation symmetries. All of these translation operations construct an Abelian translation groupRotman (1995)T. The group elements are the translation operations ^Tl1l2. For convenience, in practical calculations, the infinite system is usually replaced by a finite system with periodic boundary conditions in the in-plane directions. The dimension of the finite graphene is N1→a1×N2→a2, so integers (l1, l2) have finite possible values, i.e., l1=0,1,...,N1−1 and l2=0,1,...,N2−1.

There are N1N2 elements in the translation group T, with ^T1,0 and ^T0,1 as the two generators for this cyclic group. Using the character table of the group, it can be found that there are N1N2 one-dimensional irreducible representations for the group T. The irreducible representation can be labeled by the wave vector →k. The representation form of the translation operation is ^Tl1l2=ei→k⋅→Rl1l2 in the irreducible representation →k. For convenience, the reciprocal space is usually introduced via the definition of the reciprocal basis vectors (→b1,→b2) as,

→ai⋅→bj=2πδij,

(1)

where i and j are 1 or 2. δij is the Kronecker delta. The wave vector →k can be written in terms of the reciprocal bases as →k=k1→b1+k2→b2. Using the cyclic property of the translation group, we have ^TN1l10=1 and ^TN20l2=1, so we get k1=j1/N1 and k2=j2/N2 with j1= 0, 1, 2, …, N1−1 and j2= 0, 1, 2, …, N2−1.

The Bloch’s theorem says that, in the irreducible representation →k, the displacements of atoms in the unit cell (l1, l2) are related to atoms the (0, 0) unit cell by a phase factor ei→k⋅→Rl1l2; i.e., u(l1l2s)=u(00s)ei→k⋅→Rl1l2. A general displacement can be expanded in follow terms,

→u(l1l2s)

=

1√ms1√N1N2∑→k6∑τ=1^Q(τ)→kei→k⋅→Rl1l2→ξ(τ)(→k|00s).

This formula works for the lattice using the translational part of the space group. For nanotubes with screw symmetries (line group), the Bloch’s theorem should be generalized to include screw operations.Miloevic and Damnjanovic (1993); Popov et al. (1999); Dobardzic et al. (2003); Jiang et al. (2006); Dakic et al. (2009) An explicit introduction on the lattice dynamics of nanotubes can be found in the book chapter by Tang, Wang, and Su.Tang et al. (2011) In the above expansion, →ξ(τ)(→k|00s) is actually the vibrational displacement for atom (00s) in the →k mode (irreducible representation). These eigen vectors are orthogonal to each other. The wave vector →k is used to denote the phonon modes; τ runs over the six branches. There are two relations for the eigen vector →ξ(τ)(→k|00s) and the quantum operator ^Q(τ)−→k:

→ξ(τ)(−→k|00s)=→ξ(τ)(→k|00s)∗

^Q(τ)−→k=^Q(τ)†→k.

Dynamical Matrix

The general potential energy of graphene is determined by positions of all atoms (l1,l2,s), i.e., V0=V(→r00s,...→rl1l2s,...→rN1−1N2−1s). There will be some variance in the total potential energy, if there is a small displacement of atom (l1,l2,s), i.e., →rl1l2s→→rl1l2s+→ul1l2s. Here →ul1l2s is a small displacement for atom (l1,l2,s). The total potential energy can be expanded in a Taylor series of displacements →ul1l2s,

V(→rl1l2s+→ul1l2s)

(3)

=

V0+N1∑l1=1N2∑l2=1∑s=A,B∑α=x,y,z∂V∂uαl1l2suαl1l2s

+

12N1∑l1,l′1=1N2∑l2,l′2=1∑s,s′=A,B∑α,β=x,y,z∂2V∂uαl1l2s∂uβl′1l′2s′

⋅

uαl1l2suβl′1l′2s′+...,

where higher order nonlinear terms have been omitted. The first term V0 is the minimum potential energy of the system at the optimized configuration. The second term vanishes due to equilibrium condition for the optimized structure. The third term describes a harmonic energy induced by the displacement (vibration). Coefficients in front of the third term are normally called the force constant matrix Kl1l2sα;l′1l′2s′β=∂2V∂uαl1l2s∂uβl′1l′2s′. Applying the Bloch’s theorem to all displacement vectors, we get the potential variance up to the second order,

δV

=

12∑ττ′∑→k^Q(τ)†→k^Q(τ′)→k∑ss′=A,B∑α,β=x,y,z

(4)

⋅

Dsα;s′β(→k)ξ(τ)α(→k|00s)⋆ξ(τ′)β(→k|00s′).

We have introduced the dynamical matrix,

Dsα;s′β(→k)

=

1√msms′N1∑l1=1N2∑l2=1K00sα;l1l2s′βei→k⋅→Rl1l2,

where the summation over (l1,l2) can be truncated to the summation over neighboring atoms in case of short-range interactions. The eigenvalue of the dynamical matrix gives the eigen frequency,

∑s′βDsα;s′β(→k)ξ(τ′)β(→k|00s′)

=

ω(τ)2(→k)ξ(τ′)α(→k|00s).

(6)

The dynamical matrix is Hermitian, since the force constant matrix K is symmetric. As a result, all eigen values ω2 are real. However, there is no guarantee for the positive definiteness of the dynamical matrix, i.e., it is possible to encounter ω2<0. This positive definite property can actually be used to analyze the structure stability. The relation of ω2<0 leads to an imaginary frequency of the phonon mode, i.e., ω=iγ with real number γ. If atoms in the system are displaced according to the morphology of this imaginary mode, then the oscillation amplitude is proportional to e−iωt=eγt⟶+∞ in the limit of infinite time. An infinite vibration amplitude indicates the instability of the configuration. As an example, we note that this technique has been applied to predict the instability of nanowires,Peelaers et al. (2009) the tension induced instability of graphene,Liu et al. (2007) or the compression-induced buckling of the single-layer molybdenum disulphide.Jiang (2014a)

From the above, it is clear that the frequency from the dynamical matrix is a linear property, because it is extracted from the harmonic term. Phonon modes have infinite life time in the linear regime. The phonon life time can be limited by various scattering mechanisms. For instance, the phonon-phonon scattering becomes more important at high temperature, leading to a frequency shift and a finite value for the phonon life time. This nonlinear information can be accounted through mode coupling theoryKubo et al. (1992); Wang and Li (2004), effective phonon conception,Bonini et al. (2007); Mariani and Von Oppen (2008) or Boltzmann equation description.Ziman (1960)

ii.3 Origin of Flexural Mode

Valence Force Field Model

Figure 1: Graphene structure. (a) Honeycomb lattice of graphene. Bases →a1 and →a2 are displayed by two short blue arrows. The long red arrow illustrates a lattice vector →R−6,15=−6→a1+15→a1. (b) Sketch of the local environment of atom 1; i.e., three first-nearest and six second-nearest neighboring carbon atoms.

There are two major ingredients in the dynamical matrix. The first one is a phase factor, which is contributed by the space group symmetry of the system. As we known, all unit cells can be repeated by the (0,0) unit cell via a corresponding symmetric operation from the space group. The phase factor, ei→k⋅→Rl1l2, carries the relationship between the vibration displacement of the (0,0) unit cell and the (l1,l2) unit cell.

The second ingredient in dynamical matrix is the force constant matrix K00sα;l1l2s′β. The force constant matrix can be obtained mainly through three approaches; i.e., first-principles calculations, empirical potential, or the force constant model. For ionic materials, the shell modelCochran (1959) or the bond charge modelMartin (1969) can be useful for the description of the charge interaction.

In the first two methods, the force constant matrix is calculated by Kl1l2sα;l′1l′2s′β=∂2V∂uαl1l2s∂uβl′1l′2s′, where V is the total potential energy from an empirical potential or the Columb interaction in the first-principles calculations. uαl1l2s is the displacement of the degree of freedom (l1l2sα). This formula is realized numerically by calculating the energy change after displacing a small value for the degrees of freedom (l1l2sα) and (l′1l′2s′β). There are some existing packages for such numerical calculation. For instance, some common empirical potentials have been implemented in the lattice dynamic properties package GULP.Gale (1997) It gives the force constant matrix or the phonon dispersion directly. The first-principles package SIESTASoler et al. (2002) also gains some success in the calculation of the force constant matrix or phonon dispersion.

For the third method, there are two popular force constant models for the force constant matrix in graphene, i.e., mass-spring modelJishi et al. (1993) and valence force field model (VFFM).Aizawa et al. (1990) In the mass-spring model, each atom is denoted by a mass, that is connected to other atoms via springs. The force constant of the spring governs the force constant matrix. This model includes two-body interaction. It has been shown that the fourth-nearest neighbors should be included in the calculation of the phonon dispersion of graphene.Saito et al. (1998)

The VFFM aims to capture contribution from the valence electrons on the vibration frequency. Energy variations corresponding to both bond length and bond angle are included in this model. As pointed out by Yu in his book,Yu (2010) a big advantage of the VFFM is its transferability; i.e., force constant parameters in the model are almost the same for the same bonds within different materials. We will illustrate the explicit form of a VFFM for graphene in the following. The VFFM was successful in diamondMusgrave and Pople (1962) and CdS.Nusimovici and Birman (1967) A simplified version has been proposed by Keating,Keating (1966) which has gained success in many covalent semiconductors.

There are five VFFM terms corresponding to five typical vibration motions for the graphene sheet.Aizawa et al. (1990) The equilibrium position for atom i is →ri. The vector pointing from atoms i to j is →rij=→rj−→ri. The distance between atoms i and j is the modulus rij. We will write out interaction for one bond or one angle explicitly. Fig. 1 (b) illustrates the three first-nearest-neighboring atoms (2-4) and six second-nearest-neighboring atoms (5-10) for atom 1. The interaction for other bonds or other angles can be obtained analogously. The general expressions can be found in Refs. Aizawa et al., 1990; Jiang et al., 2006, 2008.

(1) The bond stretching interaction between atoms 1-2,

Vl=kl2[(→u2−→u1)⋅→el12]2.

(7)

kl is the force constant parameter. →el12=→r12/|→r12| is a unit vector from atom 1 to atom 2. This is the bond stretching interaction between two first-nearest-neighboring atoms. There are similar interactions for other first-nearest-neighboring carbon-carbon bonds, i.e., bond 1-3, 1-4, 2-5, and 2-6.

(2) The bond stretching interaction between atoms 1-5,

Vsl=ksl2[(→u5−→u1)⋅→el15]2

(8)

with ksl the corresponding force constant parameter. This term describes the bond stretching interaction between two second-nearest-neighboring atoms. There are similar interactions for other second-nearest-neighboring carbon-carbon bonds.

(3) The angle bending interaction for ∠213 is VBB,

VBB=kBB2(cosθ′213−cosθ213)2;.

(9)

where θ213 is the equilibrium angle and θ′213 is the angle in vibration. This interaction term describes the bending of angles, which are formed by two first-nearest-neighboring C-C bonds. There are similar interactions for the other angles: ∠213, ∠214, ∠314, ∠125, ∠126, and ∠526.

(4) The out-of-plane bond bending is a four-body interaction. It describes the interaction between atom 1 and its neighboring atoms 2-4. If atom 1 moves out of the plane, then its neighboring atoms 2-4 will try to drag it back to the plane. This potential is, Vrc,

Vrc=krc2[(3→u1−(→u2+→u3+→u4))⋅→ez]2.

(10)

→ez is the unit vector in the out-of-plane direction. Similar interaction is also applied to atom 2.

(5) If the carbon-carbon bond 1-2 is twisted, then the following twist potential will try to react the twisting motion,

Vtw=ktw2[(→u3−→u4−(→u6−→u5))⋅→ez]2.

(11)

There are similar twisting interactions for the other first-nearest-neighboring carbon-carbon bonds: 1-3, 1-4, 2-5, and 2-6.

Rigid Translational and Rotational Invariance

Figure 3: Vibration displacement for two flexural modes in graphene. Arrow on top of each atom represents the vibration component of the atom in this vibration mode. Red circles enclose small pieces of graphene, which are effectively rotated around y-axis.

The crystal is rigid in the sense that its total potential energy should not vary if the system is rigidly translated or rotated.Born and Huang (1954) According to this requirement, the empirical potential energy should satisfy two conditions.

The rigid translational invariance. It says that, if →ui=→u0 is a constant vector for all atoms, then we should have δV=0.

The rigid rotational invariance. It says that, if the system is rotated by →ui=δ→ω×→ri, then we should also have δV=0. Here, the rotation angle is |δ→ω| and the rotation direction is δ→ω|δ→ω|.

We can check that the five terms in the above VFFM satisfy both translational and rotational invariance.

For the translational invariance, the following relationship can be easily found,

→ui

=

→uj

(12)

→ui−→uj

=

0.

(13)

As a result, the translational invariance is satisfied, i.e.,

Vl=Vsl=VBB=Vrc=Vtw=0.

(14)

We will now illustrate the rigid rotational invariance for the above five VFFM potential terms.Jiang et al. (2006) During a rigid rotation motion, the displacement for atom i is

As a result, we find that (7) and (8) are zero under a rigid rotational motion, i.e.,

Vl=Vsl=0.

(18)

For the the other three potential terms (9), (10) and (11), they become summation over the following expressions, if the system is rotated rigidly,

VBB

∼kBB4[δ→ω⋅(→el12×→el13+→el13×→el12)]2=0;

(19)

Vrc

∼krc2[δ→ω×(→r12+→r13+→r14)⋅→ez1]2=0;

(20)

Vtw

∼ktw2[δ→ω×(→r43−→r56)⋅→ez12]2=0.

(21)

Fig. 2 shows the phonon dispersion in graphene along high symmetric Brillouin line ΓKMΓ using three models. The dispersion in panel (a) is calculated from the Brenner potential.Brenner et al. (2002) In particular, the lowest branch around Γ point is the flexural branch with a parabolic spectrum. Panel (b) shows that the spectrum of flexural mode is also parabolic using the VFFM. Parameters in the VFFM are kl=5.8337 eVÅ−2, ksl=5.2936 eVÅ−2, kBB=10.2245 eV, krc=14.8 N m−1, and ktw=6.24 N m−1. However, panel (c) shows that phonon spectrum of the flexural mode from the mass spring model is linear instead of parabolic. It is because the rigid rotation symmetry is violated in the mass spring model. The longitudinal and transverse parameters are considered up to the fourth nearest-neighboring atoms.Jishi et al. (1993) These parameters for the mass spring model are (kl,k⊥)= (27.7521, 4.4753) eVÅ−2, (6.8350, 0.1728) eVÅ−2, (0.5054, 0.1021) eVÅ−2, and (0.2665, 0.0012) eVÅ−2. The comparison in Fig. 2 demonstrates that the parabolic spectrum for the flexural mode is closely related to the rigid rotational invariance.

The vibration displacement of the flexural mode is actually directly related to the rigid rotational motion style. Fig. 3 shows the vibration displacement of two flexural modes. The arrow on top of each atom represents the vibration displacement of each atom in this mode. We can divide the system into lots of small pieces along the x-axis. It can be shown that each piece is effectively rotated around the y-axis in the flexural mode. Red circles in the figure illustrate two graphene pieces, which are effectively rotated around the y-axis. The rigid rotational invariance leads to zero recovery force for this vibration to first order. That is the micro-origin for the parabolic dispersion of the flexural mode.

For nanotubes with cylindrical hollow structure, the rigid rotational invariance guarantees both the zero frequency of twisting mode and the existence of flexural mode.Mahan and Jeon (2004) The flexural mode (parabolic dispersion) turns into acoustic mode (linear dispersion) gradually with the increase of the thickness of the thin plate. We have used few-layer graphene as an example to show this dimensional crossover phenomenon.Jiang and Wang (2010a)

All three of the phonon spectra in Fig. 2 are computed based on short-range empirical interaction potentials, which are known to be less accurate than first-principles calculations.Dubay and Kresse (2003); Mounet and Marzari (2005); Gillen et al. (2009) This may cause some errors because the long-range interactions may impact the flexural mode. For instance, the long-range dipolar interaction force was considered in a model for the flexural mode in graphene.Metlov (2010) In particular, the crossover in the two in-plane high frequency optical dispersion deviate from the first-principles results.

Although the phonon spectra from empirical potentials less accurate than first principles methods, they have been widely used in practice for many physical phenomena.Al-Jishi and Dresselhaus (1982); Mohr et al. (2007) For instance, classical molecular dynamics simulations are commonly used for the study of the thermal transport in graphene. In principle, the interatomic force can be calculated from first-principles calculations, but the associated computational cost renders such approaches infeasible for systems larger than a thousand or so atoms. Hence, the interatomic force is usually computed from an efficient empirical potentials like the Brenner potentialBrenner et al. (2002) or the Stillinger-Weber potential.Stillinger and Weber (1985) The thermal conductivity obtained from the molecular dynamics simulations can be explained by the phonon spectra of graphene, which should also be calculated based on the same empirical potentials for consistency. In this sense, the phonon spectra that are calculated from empirical potentials are certainly useful, though they are not as accurate as from first-principles calculations.

The inaccuracy in the optical branches in the phonon spectra from empirical potentials will impact some of the computed properties of graphene. For instance, the thermal conductivity of pure graphene is mainly limited by phonon-phonon scattering at room temperature. A typical phonon-phonon scattering process requires the involvement of the optical phonon, so the inaccuracy in the optical branches will lead to some influence on such scattering processes. However, these empirical potentials can give accurate acoustic branches in the phonon spectra, as can be seen from Fig. 2. Acoustic phonon branches are important for many physical phenomena, such as thermal transport. The thermal conductivity in graphene is mainly contributed by its three acoustic phonon branches. As a result, the inaccuracy in the optical branches in the phonon spectra has a much smaller effect on the computed values for the thermal conductivity.

iii.1 Introduction

Thermal transport occurs in the presence of temperature gradient. In metals, both electrons and phonons are important thermal energy carriers to deliver thermal energy. In insulators or semiconductors, phonons carry most of the thermal energy, while electrons only make limited contribution. The thermal conductivity contributed by phonons is called lattice thermal conductivity. Graphene is a well-known semiconductor with zero electronic band gap. Experiments show that the electronic thermal conductivity is aroundYi˘gen et al. (2013) 10 W/m⋅K, which is less than 1% of the overall thermal conductivity in graphene.Saito et al. (2007) As a result, the electronic contribution can be safely ignored in the study of the thermal conductivity in graphene. In the following, we focus on the lattice thermal conductivity in graphene.

The thermal conductivity (κ) and thermal conductance (σ) are two related concepts that are useful in different thermal transport conditions. They are related to each other from their definitions: κ/L=σ/s, where L and s are the length and the cross-sectional area, respectively. It should be noted that for quasi-two-dimensional materials like graphene, the concept of cross-sectional area is not a well defined concept, because it is only one atom thick. It is crucial to use the same thickness in the comparison of thermal conductivity from different measurements or calculations. The thermal conductance is useful in the ballistic thermal transport regime, which is typically the primary transport mechanism in nanoscale structures. The ballistic transport also happens at very low temperature, where the phonon density is too weak for phonon-phonon scattering. During ballistic transport, each phonon mode delivers a quanta of thermal energy ℏω across the system without scattering. As a result, the thermal conductance is quantized in the ballistic regime.Schwab et al. (2000) The ballistic thermal conductance does not depend on the length of the system, because of the infinite phonon mean free path.

The diffusive thermal transport happens in large systems and/or at high temperatures. The thermal transport ability is mainly limited by phonon related scattering mechanisms in the diffusive regime. In this regime, phonon modes have finite life time and finite mean free path. In other words, the thermal energy carried by a phonon mode get dissipated with increasing distance. In the diffusive regime, the thermal conductivity is a constant with respect to the length of the structure.

Experiments showed that the thermal conductivity in few-layer graphene decreases exponentially with increasing layer number and eventually crosses over to that exhibited by bulk, three-dimensional graphite value.Ghosh et al. (2010) This dimensional crossover phenomenon has received intensive theoretical effort as intrigued by the experimental work by Ghosh et al..Singh et al. (2011); Lindsay et al. (2011); Kong et al. (2009); Zhang and Zhang (2011); Zhong et al. (2011a, b); Rajabpour and Allaei (2012); Cao et al. (2012); Sun et al. (2013); Jiang (2014b) For the single-layer graphene, all of these theoretical works have shown that the single-layer graphene has the highest thermal conductivity among all few-layer graphene systems. For few-layer graphene with layer number above two, most theoretical calculations demonstrated a monotonic decrease of the thermal conductivity with increasing layer number; while the thermal conductivity was found to be independent of the layer number in a recent molecular dynamics simulations.Jiang (2014b)

It is important for nano devices to spread the thermal energy generated during its operation. Graphene has a superior thermal conductivity, which can be very helpful in delivering heat from these devices. For instance, experiments have shown that the graphene based composites have a much higher thermal conductivity owning to graphene’s superior thermal transport ability.Yan et al. (2012); Goyal and Balandin (2012); Shahil and Balandin (2012)

In the following, we focus on the role of the flexural modes on the thermal conductivity in graphene. There have been a bunch of reviews focusing on different aspects of the thermal conductivity of graphene. Wang et al. discuss the non-equilibrium Green’s function (NEGF) approach in the calculation of the thermal conductivity for nanomaterials.Wang et al. (2008, 2013) The comparison between the thermal conductivity in graphene and other carbon materials was summarized in Ref. Balandin, 2011. Different theoretical approaches to calculating the thermal conductivity in graphene is outlined in Ref. Nika and Balandin, 2012. Several reviews have been devoted to the discussion of the anomalous thermal transport in low dimensional nanoscale systems including graphene.Cahill et al. (2003); Dhar (2008); Liu et al. (2012); Li et al. (2012); Yang et al. (2012a); Luo and Chen (2013); Marconnet et al. (2013) Some basic issues on the heat transport in microscopic level are surveyed in Ref. Dubi and Ventra, 2011 by Dubi and Ventra. Zhang and Li discuss the isotopic doping effect on thermal properties on nanomaterials, including the isotopic doping effect on the thermal conductivity in graphene.Zhang and Li (2010) Heat dissipation in nanoscale electronic device is outlined by Pop in Ref. Pop, 2010.

iii.2 Simulation of Thermal Transport

There are several available approaches for the calculation of thermal conductivity. The ballistic thermal conductance can be predicted rigorously by the NEGF approach.Ozpineci and Ciraci (2001); Mingo and Yang (2003); Yamamoto and Watanabe (2006); Wang et al. (2008) For diffusive transport, several approaches are useful, such as the NEGF method,Mingo (2006); Wang et al. (2008) the mode coupling theory,Wang and Li (2004) the Boltzmann transport equation method,Spohn (2006); Lindsay et al. (2011) the classical molecular dynamic (MD) simulation,Lepri et al. (2003) and the quantum non-equilibrium MD simulation.Wang (2007); Wang et al. (2009); Ceriotti et al. (2009a, b); Dammak et al. (2009) The thermal conductivity value can be obtained directly from some open source simulation packages, such as LAMMPS.Lammps (2012) In this section, we illustrate the simulation details for thermal transport using direct MD simulation. The thermal conductivity is determined by the Fourier law.

Fourier’s Law

Fourier’s law is a linear empirical law based on observation. It states that the heat flows from high-temperature region to low-temperature region, and that the thermal current density is proportional to the temperature gradient. We consider the in-plane heat transport along x direction in Fig. 6. In this situation, the Fourier law says,

J=−κdTdx,

(22)

where J is the thermal current divided by the cross-sectional area and κ is the thermal conductivity. The minus sign on the right-hand side indicates that the heat flows from the high-temperature region to the low-temperature region.

It should be noted that the Fourier law is valid only when the local thermal equilibrium is achieved. This linear law requires the thermal conductivity to be a constant with respect to the structure dimension. However, it has been found that the Fourier law is violated in nanomaterials. More explicitly, the thermal conductivity calculated from the Fourier law is size-dependent.Chang et al. (2008); Dubi and Ventra (2009); Yang et al. (2010); Xiong et al. (2010b); Chen et al. (2013b)

According to the Fourier law, the thermal conductivity can be extracted based on the knowledge of the thermal current density and the temperature gradient. In the following, we show how to calculate these two quantities from direct MD simulations.

Figure 4: (Color online) The thermal current for different relaxation time (τ) Top: thermal fluctuation. Bottom: the averaged thermal current.

Equations of Motion

In the classical direct MD simulation, typically, structure is divided into three regions, i.e., the high temperature-controlled region, low temperature-controlled region, and the free central region. The graphene nanoribbon shown in Fig. 6 has armchair edges. The dimension is 185×12.3 Å. The thickness of the graphene is taken to be 3.35 Å. This is the inter-layer distance in the three-dimensional graphite.Saito et al. (1998) The fixed boundary is applied in the x-direction, i.e., both left and right ends are fixed during the simulation. Periodic boundary condition is applied in the y-direction. Free boundary is applied in the out-of-plane direction.

Figure 6: (Color online) The graphene is divided into three different regions, including the two fixed ends, the two temperature-controlled regions, and the free region.

In the central region, the degrees of freedom for atom i (→ri, →vi) are controlled by the following equation of motion,

d→ridt

=

→vi,

(23)

d→vidt

=

−1mi∂V∂→ri.

(24)

The Brenner potential is used to describe the interatomic interaction in this calculation.Brenner et al. (2002) The time evolution for each atom can be obtained by solving Eq. (24) numerically.

In the temperature-controlled regions, the motion of the atom is influenced by the heat bath besides the inter-atomic force. The heat bath helps to keep a constant temperature for these regions. The heat bath is described by the thermostat degrees of freedom. There are various thermostat algorithms for temperature controlling, such as Nóse-Hoover heat bath,Nose (1984); Hoover (1985) the classical or quantum Langevin heat bath,Wang (2007); Wang et al. (2009); Ceriotti et al. (2009a, b); Dammak et al. (2009) amongst others. We take the Nóse-Hoover heat bath as an example in the following demonstration of the simulation for the thermal transport. The movements of these atoms in the temperature-controlled regions are governed by the following coupled dynamic equations,

d→ridt

=

→vi,

(25)

d→vidt

=

−1mi∂V∂→ri−ηi→vi,

(26)

dηidt

=

(∑jmjv2j−gkBT)/Q.

(27)

Q

=

gkBTτ2,

(28)

g is the total degrees of freedom in the temperature-controlled region. τ represents the interaction strength between the heat bath and the system, so it is a kind of thermal relaxation time for the Nóse-Hoover heat bath. This thermal relaxation time has some direct effect on the thermal current across the system as shown in Fig. 4. The top panel is for the thermal fluctuation and the bottom panel is the average thermal current. For τ=4.5 ps, the response from the heat bath is very slow, so it corresponds to a weak coupling between the heat bath and the system. In this case, a longer thermalization time is needed to realize a stable temperature distribution across the system, although the influence from the heat bath is smaller. For τ=0.4 ps, the heat bath can respond very fast, so the heat bath couples with the system strongly. As a result, shorter thermalization time is required, but the heat bath will induce more influence to the system. Similar thermal current are obtained for τ=0.4 ps and 1.4 ps. Simulations with different relaxation time result in almost the same thermal conductivity.

Thermal Current and Temperature Gradient

The fundamental effect of the heat baths is to inject thermal energy into the system through the high-temperature region, and pump out the same amount of thermal energy in the low-temperature region. Due to energy conservation, the thermal current across the system should equal to the energy exchange between the heat bath and system in the temperature-controlled regions, as long as there is no energy accumulation in the system.Poetzsch and B¨ottger (1994); Jiang et al. (2010a)

The left-hand side is the rate of change in the total energy. As a result, we get the energy flowing from the heat bath to the system,

Eex

=

−∫t0+Δtt0∑iηimiv2idt.

(31)

Δt is the total simulation time. The summation index i runs over all atoms in the temperature-controlled region. We can thus compute the thermal current in a more symmetric manner,

J=1sEhighex−Elowex2Δt,

(32)

with s as the cross-sectional area. Ehighex and Elowex are energy from heat bath on the left and right sides. The thermal current is shown in Fig. 4. The thermal fluctuation is dJ=Ehighex−Elowex2Δt. The cross-sectional area is not included for these data shown in the figure. For τ= 0.4, 1.4, and 4.5 ps, the thermal currents are 0.87, 0.86, and 0.66 eV/ps, respectively.

Eqs. (24) and (28) are solved iteratively, by discretizing the time with a small time step, where the time step in MD simulations for graphene is usually on the order of femtoseconds. The highest-frequency phonon mode in graphene is the in-plane optical mode, with frequency around 300 THz. For a time step of 1.0 fs, there are about 20 simulation steps within one oscillation cycle of the optical mode. The trajectory of each atom and the thermostat parameter from the iterative solution are used to compute the thermal current.

The average kinetic energy for each atom gives the temperature for the atom according to the equipartition theorem 12(v2x+v2y+v2z)=32kBT, where kB is the Boltzmann constant. Hence, the temperature profile T(→ri) is obtained simultaneously from the MD simulation. Fig. 5 shows a typical temperature profile in graphene. The temperature profile within the central region x∈[0.25L,0.75L] is linearly fitted to give the temperature gradient within the system. The resulted temperature gradients are -0.28, -0.26, and -0.22 K/Å, for the three different relaxation times τ= 0.4, 1.4, and 4.5 ps. The thermal conductivity can then be extracted through the Fourier law in Eq. (22). The obtained thermal conductivity values are 120.9, 128.7, and 116.7 W/m⋅K.

iii.3 Contribution from Flexural Mode

Long Lifetime for Flexural Mode

It was shown in several recent works that the Fourier law is not valid in nanoscale structures,Cahill et al. (2003); Dhar (2008); Liu et al. (2012); Xu et al. (2014) where the thermal conductivity becomes size-dependent and increases with increasing length. For quasi-one-dimensional nanostructures, the thermal conductivity can be written as a power function of the length, i.e., κ∝Lβ, where the exponent β=0 for purely diffusive transport and β=1 for purely ballistic transport. The exponent deviates from 0, and becomes size-dependent in nanomaterials; i.e., the Fourier law is violated.

The contribution from each phonon mode to the thermal conductivity can be collected as follows,Hone et al. (1999); Gu and Chen (2007); Jiang and Wang (2011b)

κph=1V∑→kτσ→kCph(ω)v2→k.

(33)

V is the volume of the system. Cph=kBx2ex/(ex−1)2 is the heat capacity. x=ℏω/(kBT). v→k is phonon group velocity in the thermal flow direction. The lifetime for each phonon in this formula can be obtained using the single mode relaxation time approximation. Following formula gives the the lifetime for phonon mode →k due to phonon-phonon scattering,

1τps

=

(43ρL)(ℏωγ2v2z)′∑n′σ′1vgω′ω′′N(ω′,ω′′).

(34)

In this formula, ρL is the mass per length. vz is the phonon velocity along the thermal current direction. vg=|v′−v′′| is the group velocity. The energy and momentum conservation is implicated by the prime over the summation.

The phonon-phonon scattering is weak at low temperature, so the boundary scattering becomes more important, especially for systems with small size. Hence, it is also important to consider the boundary scattering process,Ziman (1960)

1τbs

=

vσ→kL×1−p1+p,

(35)

where p is the spectacular parameter. The overall phonon lifetime can be obtained as,

1τtot

=

1τps+1τbs.

(36)

The above formula gives the phonon lifetime and thermal conductivity due to boundary scattering and phonon-phonon scattering. It was found that, in the frequency range [50, 80] cm−1, the lifetime for the flexural mode is dominated by the boundary scattering, as the phonon-phonon scattering is weak.Jiang and Wang (2011b) These phonons with long lifetime transport across the system almost ballistically, while the other phonons with short lifetime behavior diffusively. As a result, the overall behavior for the thermal conductivity is sandwiched between the ballistic and diffusive transport regimes, i.e., the power factor β sits in [0,1].

The thermal conductivity in graphene was found to increase with increasing size even in the μm range.Nika et al. (2009a, b) When the width of the graphene is enlarged by a factor of 3, the thermal conductivity increases by about a factor of 1.8. There is still no conventionally accepted argument for such violation of the Fourier law in graphene.

It is clear that the three acoustic phonon branches make the largest contribution to the
the thermal conductivity for graphene. However, there is still no universally accepted fact on the relative contribution from the flexural mode to the thermal conductivity in graphene. Typically, the contribution from the flexural mode depends on the temperature or defect density or the substrate for the graphene sample. For perfect graphene, using Boltzmann transport theory, Lindsay et al. found that the flexural mode dominates the thermal conductivity of the graphene.Lindsay et al. (2011) The flexural mode contributes about 70% of the thermal conductivity in graphene, while each of the other two acoustic phonon modes contribute 10%. The large contribution from the flexural mode is attributed to the large density of states of the flexural mode, and the strict symmetry selection rule imposed on the flexural mode, which leads to very long lifetime of the flexural mode. We note that the symmetry selection rules are demonstrated in the first Brillouin zone, which was combined with the phonon-phonon scattering formula Eq. (34) to give the extremely long lifetime for the flexural mode in graphene.Nika et al. (2009a)

In other works, the flexural mode has been found to make a smaller contribution than the other two acoustic (LA and in-plane TA) modes, especially at higher temperature.Nika et al. (2009a); Aksamija and Knezevic (2011); Nika and Balandin (2012); Chen and Kumar (2012) Aksamijaa and Knezevicb found that the flexural mode has about a 50% contribution at temperatures bellow 130 K for graphene with rough edges.Aksamija and Knezevic (2011) However, the flexural mode contribution decreases quickly with increasing temperature, and the contribution from the flexural mode becomes less than 20% at 400 K. Chen and Kumar found that the LA mode dominates the thermal conductivity of both isolated and supported graphene on Cu substrate.Chen and Kumar (2012)

Furthermore, in practice, the symmetry selection rule will be relaxed in experiments, where the high symmetry of the perfect graphene is broken. For instance, the symmetry will be lowered in suspended graphene samples, which is inevitably bent during measurement. In this situation, the lifetime of the flexural mode should be considerably reduced. For deformed graphene samples, it is difficult to use the Boltzmann transport theory to compute the thermal conductivity, because the symmetry selection rule no longer valid. Therefore, classical MD simulations can be used to study the thermal transport in the deformed graphene or graphene with defects.

Several theoretical works have shown that the temperature dependence will scale as T1.5 at low temperature for the thermal conductance contributed by the flexural mode in the ballistic regime.Mingo and Broido (2005); Saito et al. (2007); Jiang et al. (2009a) In a recent experiment, Xu et al. measured the thermal conductivity of the graphene in a quasi-ballistic regime and found that the temperature dependence of the thermal conductivity scales asXu et al. (2010)T1.5, which is the same as the theoretical predictions. This consistency gives one piece of evidence that the flexural mode makes an important contribution to the thermal transport in graphene in the quasi-ballistic regime. However, it should be noted that, in practice, the sample quality plays an important role on the temperature-dependence of the thermal conductivity. In particular, the temperature-dependence for the thermal conductivity in graphene is sensitive to the defect densities or the grain size of the graphene sample.Adamyan and Zavalniuk (2012); Serov et al. (2013)

Projection Operator for Flexural Mode

We provide a normal mode projection operator for analyzing the contribution from the flexural mode to the thermal conductivity in MD simulations. From the lattice dynamic properties of the flexural mode discussed in Sec.II, we can determine the eigenvector for the flexural mode. The position of atom i is determined by the vector →ri. The eigenvectors can be used to define the normal mode projection operator, Pk,

Pk=(→ξ1,→ξ2,→ξ3,...,→ξN)

(37)

where N is the total number of atoms. From MD simulations, we have the time history of the vibration displacement of each atom,

→ui(t)=→ri(t)−→r0i.

(38)

where →r0i is the initial position for atom i. Applying the normal mode projection operator, Pk, onto the vibration displacement vector will give us a scalar normal mode coordinate, Qk(t),

Qk(t)=N∑i=1→ξ∗i⋅→ui(t).

(39)

The normal mode projection technique can selectively disclose the contribution to the thermal conductivity from each phonon mode. Fig. 7 compares the normal mode coordinate of the first flexural mode, in-plane transverse acoustic mode, and longitudinal acoustic mode from the MD simulation at 300 K. It is clear that the flexural mode has much larger normal mode amplitude than the other two modes, indicating that the flexural mode is the most important vibration morphology in the graphene during the MD simulation. It shows explicitly that the major contribution is from the flexural mode to the thermal conductivity.

Tuning the Thermal Conductivity Via Flexural Mode

The graphene is very flexurable in the out-of-plane direction, so it is easier to modify the flexural mode. Hence, it will be an efficient way to manipulate the thermal conductivity through the flexural motion.

The first method is to introduce inter-layer coupling for the flexural mode; e.g., coupling different graphene monolayers to form a few-layer graphene. Indeed, it was observed in experiment that the thermal conductivity is considerably reduced by increasing layer numbers in few-layer graphene.Ghosh et al. (2010) The thermal conductivity in bilayer graphene is about 30% lower than that of the single-layer graphene. The reduction of the thermal conductivity is attributed to the enhanced phonon-phonon scattering for the flexural mode in thicker few-layer graphene, where more scattering channels are available.

The second direct method to manipulate the flexural mode is to fold the graphene as shown in Fig. 8. A fold will introduce inter-layer coupling to the flexural mode.Ouyang et al. (2011); Yang et al. (2012b) This coupling can be further increased by compression the inter-layer space in the folds. For a flat graphene, flexural mode is difficult to scatter with phonon modes from the high-frequency optical branches, because the band gap between the flexural mode and the optical branches is very large due to the parabolic spectrum of the flexural mode. However, the phonon-phonon scattering channels are considerably increased when the inter-layer space is compressed in the folds. The compression of the inter-layer space leads to considerable shifting up of the flexural phonon spectrum, thus narrowing the band gap between the flexural mode and the optical branches. As a result, the phonon-phonon scattering channels increase significantly. This results in a reduced thermal conductivity in folded graphene.

The third method to manipulate the flexural mode is to put the graphene on a substrate.Cai et al. (2010); Chen et al. (2011); Lee et al. (2011); Guo et al. (2011); Ong and Pop (2011); Huang et al. (2011); Qiu and Ruan (2012); Guo et al. (2012) The thermal conductivity can be reduced by an order of magnitude due to damping of the flexural mode by the substrate.Ong and Pop (2011) Furthermore, the increase in the coupling strength between graphene and substrate will enhance the thermal conductivity. It is because of the coupling of the flexural mode to the substrate Rayleigh waves, which results in a hybridized mode with linear spectrum and higher group velocity than the original flexural mode in graphene. Quite recently, Amorim and Guinea examined the flexural mode of graphene on different substrates, with the consideration of the dynamics of the substrate.Amorim and Guinea (2013)

iv.1 Introduction

Negative coefficient of thermal expansion (CTE) occurs in many materials, typically at low temperature. For example, bulk Si and Ge are two well-known semiconductors with negative CTE at temperatures below 100 K, which is attributed to the inter-play between the bond-stretching and bond-bending forces.Xu et al. (1991); Wei et al. (1994) Ref. Barrera et al., 2005 summarizes the negative CTE in different materials and various possible mechanisms underlying this phenomenon.

The CTE for graphene was also found to be negative. On the experimental side, in 2009, the CTE is found to be about −7×10−6 K−1 at room temperature as measured by Bao et al.Bao et al. (2009) The room temperature CTE is about −8×10−6 K−1 in the experiment by Yoon et al..Yoon et al. (2011) The measured CTE is negative in a wide temperature range. Singh et al. found that graphene has a negative CTE for temperatures bellow 300 K,Singh et al. (2010) while Yoon et al. obtained a negative CTE for graphene for temperatures bellow 400 K.Yoon et al. (2011)

On the theoretical side, there are several standard approaches to compute the CTE. The lattice constant is temperature dependent, which can be used to extract the CTE value.Mounet and Marzari (2005); Zakharchenko et al. (2009); Pozzo et al. (2011) These calculations show a minimum lattice constant at a transition temperature. Below this transition temperature, the CTE decreases with increasing temperature. Above this transition temperature, the CTE increases with increasing temperature. This implies that CTE is negative below the transition temperature and becomes positive above the transition temperature. In our recent work, the NEGF approach was implemented to compute the CTE of graphene.Jiang et al. (2009c) As an advantage of the NEGF approach, this method is able to account for the quantum zero-point vibration effect in a quite natural manner. Another advantage of the NEGF approach is that it is straightforward to decompose the contribution of each phonon mode to the total CTE value. We will mainly discuss this NEGF approach in the following.

iv.2 Green’s Function Approach for Thermal Expansion

The thermal expansion is usually studied by MD simulations or the standard Grüneisen method.Gruneisen (1926) In our recent work, the NEGF approach is used to examine the thermal expansion phenomenon in the single-walled carbon nanotube and graphene, where the calculated CTE agrees quite well with experiments.Jiang et al. (2009c) In this section, we review some key steps in this NEGF approach for the thermal expansion.

Green’s Function Approach

The potential energy associated with the phonon modes are assumed to be,

V

=

∑ijKij2uiuj+Hn.

(40)

It includes linear and nonlinear interactions. Both third and fourth orders nonlinear interactions are considered,

Hn=∑lmnklmn3ulumun+∑opqrkopqr4uoupuqur.

(41)

Klmn and Kopqr are force constant matrix elements. We have extracted the nonlinear coefficients klmn and kopqr from the Brenner potential.Brenner et al. (2002)

Details for the Green’s function approach in thermal properties can be found in Refs. Wang et al., 2008, 2013. We utilize the following two GFs for the study of the thermal expansion,

Gj(τ)

=

−iℏ⟨TτuHj(τ)⟩,

(42)

Gjk(τ,τ′)

=

−iℏ⟨TτuHj(τ)uHk(τ′)⟩.

(43)

uHj(τ) is the vibrational displacement. For convenience, uHj(τ) also includes the square root of the atom’s mass. It is very nice that Gj can be obtained analytically,

Gj=∑lmnklmnG>lm(0)~Grnj[0].

(44)

G>(0) is the greater GF in time domain. ~Gr[0] is the retarded GF in frequency domain. These two GFs can be computed without any integration. Using Eq. (42), we get the average vibrational displacement for atom j, i.e., ⟨uj⟩. CTE can be computed from the derivative of the displacement with respect to the temperature.

It is more convenience to work in the normal mode space for the derivative of the one-point GF. After the Fourier transformation, we get,

where f is the Bose distribution function. K is the force constant matrix. S stores the eigen vector of K, i.e., S†KS is diagonal with diagonal elements ω2μ.

A particular boundary condition is used for the NEGF treatment of the thermal expansion phenomenon. Specifically, the left end is fixed, while the right end is free. Periodic boundary condition is applied in the in-plane lateral direction. It was shown that graphene has a negative CTE for temperatures bellow 600 K.Jiang et al. (2009c) The minimum CTE value is achieved around 200 K. These results are in good agreement with the experimental findings. The NEGF approach is able to examine the substrate coupling effect, which has been found to be important for the CTE value of graphene.Pozzo et al. (2011)

Green’s Function Approach and Grüneisen Method

We will show that the NEGF approach is equivalent to the traditional Grüneisen method in the weak nonlinear limit, if the system is isotropically and uniformly deformed in the thermal expansion phenomenon.Wang and Jiang (2014)

From above, the CTE from NEGF for atom N is,

αN

=

1L∑lmnμklmnSlμS∗mμcμω2μ(−K−1nN),

(46)

where cμ=ℏωμdf/dT=Crmph(ωμ) is the heat capacity for phonon mode μ as introduced in Eq. (33). The Grüneisen parameter for mode μ is,

γμ

=

−∂lnωμ∂lnV

(47)

=

−L3ω2μ∑lmnSlμS∗mμklmnϵn,

where ϵn=δRn/L is the change of position for atom n with respective to total length L. The CTE from Grüneisen method is then,

α

=

13B(∂P∂T)V

(48)

=

13BV∑μγμcμ

=

−1B∑lmnμϵnklmnSlμS∗mμcμω2μ,

where P is pressure, B is bulk modulus, V is volume, and an isotropic and uniform deformation has been assumed for the system during the thermal expansion phenomenon.

To show the equivalence between these two methods, let’s assume the system undergoes a uniform deformation by external force F. In this situation, the displacement of atom n can be obtained from the force constant matrix, or equivalently from the bulk modulus,

un

=

K−1nNF=ϵnFLB,

(49)

so,

1Bϵn

=

K−1nNL.

(50)

As a result, we find that,

α

=

αN.

From the comparison, we show that the NEGF approach has three advanced properties over the Grüneisen method. First, the NEGF approach can be applied to study structures which have a nonuniform deformation, while the Grüneisen method can only be applied for structures with uniform deformation in the thermal expansion phenomenon. Second, the NEGF approach is applicable for nanostructures without periodicity, while the Grüneisen method requires periodicity with a well-defined bulk modulus. Third, the NEGF approach is suitable for systems with anisotropic CTE, while Grüneisen method only works for systems with isotropic CTE.

In the meantime, we point out two disadvantages for the NEGF approach as compared with the Grüneisen method. First, the calculation of K−1 is computationally expensive in the NEGF approach, especially for large systems. Hence, in terms of computation cost, the Grüneisen method is better than the NEGF approach for materials with isotropic CTE, where both methods are applicable. Second, the NEGF approach is based on a perturbation theorem and high-order nonlinear terms have been omitted; while Grüneisen method includes an overall effect from all high-order nonlinear interaction terms.

iv.3 Contribution from Flexural Mode

Figure 9: Contribution to CTE from six lowest-frequency phonon modes in graphene sheet without substrate interaction. Blue solid line is for graphene with (length, width)=(10, 8.5) Å, green dashed line is for (20, 17) Å and red dotted line is for (20, 8.5) Å. Insets are the vibrational morphology for the corresponding mode. (a), (b), (c) and (e) are the first four bending modes. (d) is a tearing mode. (f) is the longitudinal phonon mode.

We will now demonstrate that the large negative CTE in graphene is due to its flexural modes. Eq. (45) includes the contribution from all phonon modes to the CTE through the diagonal matrix,

⎛⎜
⎜
⎜
⎜
⎜⎝⋱1ωμ(dfdT)⋱⎞⎟
⎟
⎟
⎟
⎟⎠.

This matrix is diagonal, which indicates that each phonon mode makes separate contributes for the negative CTE. If this matrix has a single nonzero element, 1ωμ(dfdT), then Eq. (45) yields the sole contribution from the mode μ. Hence, we can distinguish the independent contribution from each phonon mode to the CTE.

Fig. 9 shows the contribution from the six lowest-frequency phonon modes. There is no substrate interaction in this figure, i.e., γ=0. The inset in each panel is the corresponding vibration morphology of the phonon mode. (a), (b), (c) and (e) are the first four bending modes and (d) is an interesting tearing mode. Among all of these six phonon modes, the first bending mode shown in (a) has 90% contribution to CTE. Due to its bending morphology, this mode induces a contraction effect in the graphene sheet, and thus the negative CTE.

If the substrate interaction is nonzero as shown in Fig. 10, the CTE is clearly enhanced and the size effect becomes weaker. The contribution from the second bending mode is also very important. The substrate interaction is more important for larger piece of graphene, because the bending movement in larger graphene is more serious than smaller piece of graphene when it is fixed on one edge. In case of strong substrate interaction, graphene is so difficult to be bent that all bending modes do not contribute. In this situation, the CTE is dominated by the sixth mode, which leads to a positive CTE. As a result, the CTE is positive in whole temperature range, and the size effect on the CTE becomes smaller.

v.1 Introduction

Graphene has many remarkable mechanical properties. Readers are referred to Refs. Li et al., 2009; Terrones et al., 2010; Yang et al., 2012c; Lau et al., 2012; Tu and Ou-Yang, 2008 for comprehensive reviews on various mechanical properties of graphene. For instance, the Young’s modulus for graphene is on the order of TPa. There are several different approaches for the investigation of the Young’s modulus. The atomic force microscope can measure the force-displacement relationship for graphene, from which the Young’s modulus can be extracted. This method is used in the experiment in 2008, which found the Young’s modulus for graphene to be 1.0±0.1 TPa. This method has also been adopted in some theoretical works to study the Young’s modulus in graphene.Lu (1997); Lier et al. (2000); Kudin et al. (2001); Konstantinova et al. (2006); Reddy et al. (2006); Huang et al. (2006); Liu et al. (2007); Khare et al. (2007); Wei et al. (2009); Jiang et al. (2010d); Yi and Chang (2012); Zhou et al. (2013)

The classical theory of elasticity has been widely used to predict the mechanical properties for graphene. Finite element simulations, which are based on a numerical discretization of the equations of elasticity, were used to explain the edge stress induced warping phenomenon in the graphene nanoribbon.Shenoy et al. (2008) The finite element method has also been coupled with atomistic calculations to simulate the crack propagation in graphene in an efficient way.Khare et al. (2007) The elasticity approach has also been combined with atomic potentials to study elastic properties for finite size graphene.Reddy et al. (2006); Duan and Wang (2009)

As an ultra-thin plate, graphene’s flexural mode is directly related to its Young’s modulus.Landau and Lifshitz (1995) The flexural mode has the lowest frequency, so this mode is the easiest to be excited by thermal vibration. As a result, the thermal vibration of graphene is closely related to its Young’s modulus. Using this relationship, it is possible to extract the Young’s modulus from the thermal vibrations. In 1996, this idea was implemented by Treacy et al. to measure the Young’s modulus of carbon nanotubes.Treacy et al. (1996); Krishnan et al. (1998) In our recent work, we applied this method to compute the Young’s modulus from graphene’s thermal vibration.Jiang et al. (2009d) In this section, we will review this approach and emphasize the relationship between the flexural mode and the Young’s modulus of graphene.

v.2 Flexural Mode and Young’s Modulus

In the elasticity theory, the flexural mode is directly related to the Young’s modulus of a plate. Here we review some key derivation steps for the Young’s modulus of graphene.Jiang et al. (2009d) The flexure mode for the elastic plate is governed by the following equation,Landau and Lifshitz (1995)

ρ∂2z∂t2+DhΔ2z=0,

(51)

where ρ is the mass density. D=112Yh3/(1−μ2) is the bending modulus. Y is the Young’s moduls. μ is the Poisson ratio. h is the thickness of the plate. This partial differential equation can be solved after the following boundary conditions are applied,

z(t,x=0,y)

=

0,

z(t,x=L,y)

=

0,

(52)

z(t,x,y+L)

=

z(t,x,y).

According to this boundary condition, we can get following solution,Polyanin (2002)