InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Nanotechnology and Nanomaterials » "Graphene - New Trends and Developments", book edited by Farzad Ebrahimi, ISBN 978-953-51-2220-3, Published: November 25, 2015 under CC BY 3.0 license. © The Author(s).

Chapter 1

A Review on Modeling, Synthesis, and Properties of Graphene

By Farzad Ebrahimi and Ebrahim Heidari
DOI: 10.5772/61564

Article top


Mechanical testing of graphene. Schematic of nanoindentation on suspended graphene membrane (left figure). Atomic force microscope image of a fractured graphene membrane (right figure).
Figure 1. Mechanical testing of graphene. Schematic of nanoindentation on suspended graphene membrane (left figure). Atomic force microscope image of a fractured graphene membrane (right figure).
Schematic diagram of the setup for the tearing studies of graphene: side and top views. The inset shows the sheet orientation. An initial flap of 8 nm in width is cut in the sheet, folded back, and moved at a constant speed.
Figure 2. Schematic diagram of the setup for the tearing studies of graphene: side and top views. The inset shows the sheet orientation. An initial flap of 8 nm in width is cut in the sheet, folded back, and moved at a constant speed.
(a) Thermal conductivity κ as a function of temperature for representative data of suspended graphene [55], SiO2-supported graphene [64], ~20-nm-wide graphene nanoribbons (GNRs)[63], suspended single-walled CNTs (SWCNTs)[66], multi-walled CNTs (MWCNTs)[67], type IIa diamond, graphite in-plane and out-of-plane. Additional data for graphene and related materials are also summarized in Ref.[54]. (b) Room temperature ranges of thermal conductivity data κ for diamond [57], graphite (in-plane) [54], carbon nanotubes (CNTs) [54], suspended graphene [54, 55], SiO2-supported graphene [64], SiO2-encased graphene [65], and GNRs [63].
Figure 3. (a) Thermal conductivity κ as a function of temperature for representative data of suspended graphene [55], SiO2-supported graphene [64], ~20-nm-wide graphene nanoribbons (GNRs)[63], suspended single-walled CNTs (SWCNTs)[66], multi-walled CNTs (MWCNTs)[67], type IIa diamond, graphite in-plane and out-of-plane. Additional data for graphene and related materials are also summarized in Ref.[54]. (b) Room temperature ranges of thermal conductivity data κ for diamond [57], graphite (in-plane) [54], carbon nanotubes (CNTs) [54], suspended graphene [54, 55], SiO2-supported graphene [64], SiO2-encased graphene [65], and GNRs [63].
(a) The resistivity of a single layer of graphene vs. gate voltage. (b) The quantum Hall effect in single-layer graphene. Figures taken from (Novoselov, et al. (2005).
Figure 4. (a) The resistivity of a single layer of graphene vs. gate voltage. (b) The quantum Hall effect in single-layer graphene. Figures taken from (Novoselov, et al. (2005).

A Review on Modeling, Synthesis, and Properties of Graphene

Farzad Ebrahimi1 and Ebrahim Heidari1

1. Introduction

Graphene is the name given to a two-dimensional sheet of sp2-hybridized carbon. Its extended honeycomb network is the basic building block of other important allotropes; it can be stacked to form 3D graphite, rolled to form 1D nanotubes, and wrapped to form 0D fullerenes. Graphene is the name given to a two-dimensional sheet of sp2-hybridized carbon [1]. Graphene is a wonder material with many superlatives to its name. It is the thinnest known material in the universe and the strongest ever measured. Its charge carriers’ exhibit giant intrinsic mobility, have zero effective mass, and can travel for micrometers without scattering at room temperature [2]. Synthesis and characterization of graphenes pose challenges, but there has been considerable progress in the last year or so [3]. A great deal of research has been conducted to explore the promising properties of the single-layered graphene sheets (SLGSs) after the appearance of the new method of graphene sheet preparation. Stankovich et al. have proposed their findings related to the synthesis and exfoliation of isocyanate-treated graphene oxide nanoplatelets [4]. Implementing the chemical reduction, they have also been able to produce the graphene-based nanosheets [5]. In addition, Ferrari has reported the Raman spectroscopy of the SLGS [6]. Furthermore, Katsnelson and Novoselov have explored the unique electronic properties of the SLGSs [7]. They have stated that the graphene sheet is an unexpected bridge between condensed matter physics and quantum electrodynamics. On the other hand, Bunch et al. have reported the experimental results of using electromechanical resonators made from suspended single- and multi-layered graphene sheets [8]. The superior mechanical, chemical, and electronic properties of nanostructures make them favorable for nano engineering applications [9]. Graphene sheets are one of the most important nano-sized structural elements which are commonly used as components in micro-electro-mechanical systems (MEMS) and nano-electro-mechanical systems (NEMS) [8, 10]. Furthermore, it has been revealed that adding graphene sheets to polymer matrix could greatly improve the mechanical properties of the host polymer [11]. In addition, nanostructures such as armchair carbon nanotubes and nanoplates have shown significant potential applications in the field of environmental technologies [12]. Nano-mechanical resonators are one of most important NEMS devices which have received increasing attention from the scientific community in recent years [13-16]. The nano mechanical resonators may operate at very high frequencies up to gigahertz range [17]. The potential applications of the SLGSs as mass sensors and atomistic dust detectors have further been investigated [10]. Also, the promising usage of the SLGS as strain sensor has been examined [18].

Since graphene has a prominent application in human’s life, the necessity of mechanic analytical approach for graphene is drastically felt. There are many approaches to analyze a graphene and other nonoplates mechanically, however, they can be divided into two bunches: first, the methods that consider graphene or other nonoplates downright, and second, the methods that consider interactions between graphene and other nonoplates with their surrounding environment.

The first cluster of approaches is often the Molecular Dynamics (MD) method. It is very powerful method has furthered scientists in their case studies. Sheehan et al. utilized the molecular dynamics methodology for analyzing the effect of solvents on reaction kinetics and post reaction separation is presented [19]. Kresse et al. used ab initio molecular dynamics to predict the wave functions for new ionic positions using sub-space alignment [20].With the ability to examine atomic-scale dynamics in great detail, researchers have used MD to gain new insight into problems that have been resistant to theoretical solution, such as solid fracture [21], surface friction [22], and plasticity [23]. For example, a 10-nm cubic domain of a metal can be simulated only for times less than around 10-10 s, even on very large parallel machines [22]. Increases in this simulation time require a proportional reduction in the number of atoms simulated. The results of such a simulation therefore can rarely be compared directly to experiments, since laboratory observations of these sorts of mechanical phenomena are usually made on much larger length and time scales. One possible approach that can be applied to many problems is to use MD only in localized regions in which the atomic-scale dynamics are important and a continuum simulation method (such as finite elements) everywhere else. This general approach has been taken by several different groups of researchers. Abraham et al. [24, 25] have developed a coupled finite element, MD, tight-binding (FE/MD/TB) method in which the three methods are used concurrently in different regions of the computational domain. Another method developed recently is the quasi continuum method [26-31], in which atomic degrees of freedom are selectively removed from the problem by interpolating from a subset of representative atoms, similar to finite element interpolation in which atomic degrees of freedom are selectively removed from the problem by interpolating from a subset of representative atoms, similar to finite element interpolation. Finally Wagner et al. [32] submitted a professional coupling of atomistic and continuum simulation method that is called bridging-scale method. In this review, we are going to represent some properties of graphene and study briefly the mechanic analytical approaches that we mentioned before.

2. Structure, synthesis, and properties

Graphene has a honeycomb network that could have ripple in the surface. Ripples can induce the local electrical and optical properties of graphene. Three different types of graphenes can be defined: single-layer graphene (SG), bilayer graphene (BG), and few-layer graphene (FG, number of layers).

Typically, the important properties of graphene are a quantum Hall effect at room temperature, an ambipolar electric field effect along with ballistic conduction of charge carriers, tunable band gap, and high elasticity. Although graphene is expected to be perfectly flat, ripples occur because of thermal fluctuations. Ideally, graphene is a single-layer material, but graphene samples with two or more layers are being investigated with equal interest.

There are now four primary ways to produce ‘pristine’ graphene:

  1. Epitaxial graphene: This method involves chemical vapor deposition (CVD growth) on epitaxially matched metal surfaces

  2. Micromechanical Exfoliation or micromechanical cleavage: in which highly oriented pyrolitic graphite (HOPG) is pealed using Scotch-tape and deposited on to a silicon substrate.

  3. Exfoliation of graphite in solvents: Gram quantities of single-layer graphene have been prepared by employing a solvothermal procedure and subsequent by sonication. In this process, the solvothermal product of sodium and ethanol is subjected to low-temperature flash pyrolysis yielding a fused array of graphene sheets, which are dispersed by mild sonication. A single-layer graphene can be produced in good yields by the solution-phase exfoliation of graphite in an organic solvent, such as N-methylpyrrolidone (NMP). This process works because the energy required to exfoliate graphene is balanced by the solvent–graphene interaction.

  4. “Other methods,”: such as

    • Substrate-free gas-phase synthesis of graphene platelets in a microwave plasma reactor

    • Arc discharge synthesis of multi-layered graphene

    • Graphene can be grown on metal surfaces by surface segregation of carbon or by decomposition of hydrocarbons.

    • etc.

2.1. Mechanical properties

Pristine graphene structures are found in 2D plane sheets. It has a hexagonal crystal lattice which resulted in covalent bonds between carbon atoms. In the environment, graphene are discovered in tow forms: “monolayer” and “free- standing.” With the first form, we find graphene parts as a cover over a substrate material such as SiC. However, we are able to find graphene individually and independent from other materials in the environment which corresponds with the second form, “free-standing graphene” [33].

Mechanical properties for any crystal material are affected by pristine Lattice and defects are comprised of dislocations and grain boundaries [34, 35]. For example, we can mention elastic properties of materials that are affected by atoms interactions and lattice geometry, whereas strength and plastic flow stress as another properties of materials are affected by characteristics of defects. Indeed caused defects in the material severely decrease strength of it in comparison of ideal material. Anyway, we are always not able to impede the existence of defect and its effect in the materials. However, there is one exception; nano-materials can be discovered defect-free initially, and this is the main reason of superiority of strength for these materials [36]. Graphene as a nano-material is not excepted in this issue.

Lee and his co-workers performed the pioneer empirical analysis of elastic properties and strength of pristine graphene [37]. A deposited graphene membrane onto a substrate material that possesses some cavities on the surface is loaded by the tip of the atomic force microscope (Fig. 1) [37], and it was discovered that graphene brings out both nonlinear elastic behavior and brittle fracture. Thus, for nonlinear elastic behavior, we can write:  σ=Eε+Dε2, where σ is the applied stress (the symmetric second Piola Kirchhoff stress), E is the Young modulus, ε is the elastic strain (uniaxial Lagrangian strain), and D is the third-order elastic stiffness. This experiment is convoyed by this result as follows: Young modulus of E = 1.0 TPa, and the third-order elastic stiffness of D = -2.0 TPa. The Young modulus they found is very close to the Young modulus of nanotubes. They also found that brittle fracture happens at an intrinsic stress as much as σint=130GPA, which is very huge and magnitude.

Simulating by computer [38] shows E = 1.05 TPa and σint=110 GPa for Young modulus and brittle fracture of graphene, which is compatible with explorations of Lee and his co-workers. All of these explorations prove that graphene can be very useful for structural applications and for the cases that we are dependent on high strength. Furthermore, graphene is flexible and can be bent easily, which make it more desirable and attractive.

The propagation of crack in monolayer graphene has been studied empirically and also analytically (molecular dynamics) by considering Crystallographic characteristics [39]. It has been evident that the sources of cracks in monolayer graphene membranes are unavoidable, mechanically applied stresses that are exerted during their processing. Cracks or tears propagate along the sides of the hexagonal crystal lattice and defray an occasional direction change as much as 30o in vertices of hexagonal.


Reprinted with permission from C. Lee, X. Wei, J.W. Kysar, J. Hone, Science, Volume 321, 385 388, 2008. Copyright (2008) by the American Association for the Advancement of Science.

Figure 1.

Mechanical testing of graphene. Schematic of nanoindentation on suspended graphene membrane (left figure). Atomic force microscope image of a fractured graphene membrane (right figure).

It is sometimes seen that propagation can go under the TEM electron beam [40]. Kim et al. [39] used this assumption that the propagation of crack is motivated by incorporating the effects of stress concentrations at the tip of the crack and the ionization influence of electron beam. Under the simultaneous effect of these problems, atomic bonds break in the vicinities of a crack tip, and propagation takes place.

Novoselov et al. [40] performed a computer simulation on the exfoliation of graphene sheets from adhesive substrate to examine crystallographic features of cracks which happens. The idea of this simulation developed from this fact that graphene ribbons became tapered as they were produced by exfoliation process (Fig. 2). They found that tear angle is affected by adhesive strength. Their simulation showed that, with the low strength adhesive, tearing takes place in the conqueror armchair direction of the hexagonal crystal lattice of graphene, and meanwhile, occasional change in direction is observed rarely. They also concluded that any increase in adhesive strength results in more tear angel; hence, in pretty high adhesive strength a change of 60o in the direction of tearing will present.

As a material, graphene is not excluded from defects. The defects that graphene may obsess are: vacancies [41], Stone Wales defects [41], dislocations [42], and grand boundaries of GBs. Among these defects, dislocations and GBs are very common and play a prominent role in mechanical properties of graphene. For instance, dislocations can cause plastic flow in graphene, whereas GBs decrease its strength characteristics [44, 45]. Dislocations can also violate translational symmetry of graphene [46].

There is no getting around GBs in graphene specimens when it is produced in a large area; thus, their effects on mechanical properties have been always noticeable from both fundamental and applied sides. Empirical data [44] have proved that free-standing polycrystalline graphene under concentrated load has severe low failure strength when found in polycrystalline state than when produced as a single crystal. In this experiment, they used a tip of atomic force microscopy (AFM) to exert a concentrated load on a polycrystalline graphene membrane, and it was observed that the force needed to cause a tear in the graphene along GBs is an amount of 100 nN [44], this is while the force for tearing a single-crystal exfoliated graphene is not more than 1.7 mN [37].

Vargas et al. [45] did similar test using atomic force microscopy and molecular dynamics simulations to study the mechanical characteristics of a graphene with polycrystalline structure that is obtained by chemical vapor deposition. They used nanoindentation measurements and found that out-of-plane ripples effectively decrease in-plane stiffness in the mentioned graphene. They also found that GBs effectively decrease the breaking strength of the graphene. Molecular dynamics showed them voids can significantly weaken the graphene membranes. In fact, GB is a place where amorphous carbon and iron oxide nanoparticles are absorbed [45].


Reprinted with permission from D. Sen, K.S. Novoselov, P. Reis, M.J. Buehler, Small, Volume 6, 1108,2010.

Figure 2.

Schematic diagram of the setup for the tearing studies of graphene: side and top views. The inset shows the sheet orientation. An initial flap of 8 nm in width is cut in the sheet, folded back, and moved at a constant speed.

2.2. Thermal Properties

The heat flow direction in a two dimensional graphene can be divided into in-plane and out-of-plane directions. In-plane heat flow is greater than out-of-plane one and is developed due to covalent sp2 bonding between carbon atoms, while the latter is emanated from weak van der Waals coupling.

Graphene transistors and interconnectors are beneficiaries of in-plane heat flow depending on a certain channel length. Although thermal coupling with substrate materials is miserly, it is a prominent reason for the dissipation of heat flow. We can regulate the heat flow by phonon scattering, edges, or interfaces. Ultimately, the unusual thermal properties of graphene stem from its 2D nature, forming a rich playground for new discoveries of heat-flow physics and potentially leading to novel thermal management applications.

By reviewing thermal properties, with no respect to material, one that should be considered is the specific heat. This is a quantity that implies two things: first, the thermal energy that a body is capable to store, and second, the rate of cooling and heating that a body will experiment. The later can be modeled by the thermal time constant τ ≈ RCV, where τ is the thermal time constant, R is the thermal resistance for heat dissipation (the inverse of conductance, R=1/G) and V is the volume of the body. Thermal time constants can be varied from 0.1 ns for a single graphene sheet or carbon nanotube (CNT) and 10 ns for nanoscale transistors to 1 ps for the relaxation of individual phonon modes [47 -49].

The specific heat of graphene is divided into phonons (lattice vibrations) and free electrons contributions, C=Cp+Ce, Knowing that phonon contributions are dominant [50 -52]. Phonon specific heat as a dominant coefficient becomes constant at very high temperature near in-plane Debye temperature ΘD ≈2100 K, At this time, we have Cp=3NAkB ≈25 J mol –1K–1 ≈ 2.1 J g–1K–1, also known as the Dulong–Petit limit, where NA is the Avogadro’ number and kB is the Boltzmann constant. This property belongs in a classical solid behavior when six atomic degrees of motion are entirely exited and each carries 1/2 kBT energy.

In heat transport exploration, it is assumed that the thickness of a graphene monolayer is about the graphite interlayer spacing h3.35 A°'. Graphene has one of the highest in-plane thermal conductivities at room temperature, about 2000–4000 W m–1 K–1, when it was found in freely suspended samples [53, 54]. This range corresponds with values between isotopically purified samples (0.01% C instead of 1.1% natural abundance) as a right end with large grains [55] and isotopically mixed samples or those with smaller grain sizes as a left end.

Disorders with no respect of their source introduce more phonon scattering, and this results in a descendant of conductivity lower than the mentioned range [56]. Figure 3a compares the thermal conductivity of natural diamond (about 2200 W m–1 K–1) with those of other related materials at room temperature [57, 58]. Figure 3b exhibits the thermal conductivity of materials in Figure 3a with respect to lack of disorders.

Heat flow is strongly limited by weak van der Waals interaction in both of directions: cross-plane (along the c axis) and perpendicular to a graphene sheet, knowing that the van der Waals interaction in the perpendicular direction is between graphene and adjacent substrates, such as Sio2. As we can see in Figure 3a, the thermal conductivity along the c axis of pyrolytic graphite is just ~6 W m-1 K-1 at room temperature. The relevant metric for heat flow across such interfaces is the thermal conductance per unit area, G″ = Q″ / Δ T ≈ 50 MW m –2 K –1 at room temperature [60 -62]. This is approximately equivalent to the thermal resistance of an ∼ 25-nm layer of SiO2 [59] and could become a limiting dissipation bottleneck in highly scaled graphene devices and interconnects [63], as discussed in a later section. When we have a few layers of graphene (from 1 to 10 layers), it can be expertized that interlayer resistance, 1/ G″, remains almost constant and pretty smaller than the resistance between graphene and its nongraphene environment [61]. Indeed, the interlayer thermal conductance of bulk graphite is ∼ 24 GW m –2 K –1 if the typical 3.35- A° spacing and the c axis thermal conductivity are assumed.

It must be remarked that surface effects are able to decrease the thermal conductivity of graphene because of the sensitivity of phonon propagation to surface or edge perturbation, and as a result of this, the in-plane thermal conductivity of freely suspended graphene is drastically lower than a graphene nanoribbon or a graphene contacted with a substrate.


Figure 3.

(a) Thermal conductivity κ as a function of temperature for representative data of suspended graphene [55], SiO2-supported graphene [64], ~20-nm-wide graphene nanoribbons (GNRs)[63], suspended single-walled CNTs (SWCNTs)[66], multi-walled CNTs (MWCNTs)[67], type IIa diamond, graphite in-plane and out-of-plane. Additional data for graphene and related materials are also summarized in Ref.[54]. (b) Room temperature ranges of thermal conductivity data κ for diamond [57], graphite (in-plane) [54], carbon nanotubes (CNTs) [54], suspended graphene [54, 55], SiO2-supported graphene [64], SiO2-encased graphene [65], and GNRs [63].

It has been seen that the in-plane thermal conductance G of graphene can reach a significant fraction of the theoretical ballistic limit in sub-micrometer samples, owing to the large phonon mean free path (λ ≈ 100 to 600 nm in supported and suspended samples, respectively). However, thermal properties of graphene could be highly tunable, so that makes it useful for heat- sinking applications when we regulate it in ultra-high thermal conductivity, and it is useful for thermoelectric applications when it is regulated for ultra-low thermal conductivity.

2.3. Electrical Properties of Graphene

The electronic properties of graphene are one of the prominent cases that relating experimental researches have dealt with. The controllable continuous transformation of charge carriers from holes to electrons was one of the most notable features in pioneering researches.

An example of the gate (or gate electrode is a region at the top of the transistor whose electrical state determines whether the transistor is on or off) dependence in single-layer graphene is shown in Fig. 4a. This dependency is much weaker in multiple layers of graphene because electric field is screened by the other layers.

The high electronic mobility of graphene permits the development of quantum hall effect (an effect that is visible in conductor and semiconductors materials; when there are both electrical and magnetic field at the same time in a conductor or semiconductor material, it can arise an electric potential perpendicular to the magnetic field that causes electric current perpendicular to the magnetic field) at low temperatures and high magnetic fields, for both electrons and holes (Fig. 4b) (Novoselov, et al., 2005; Zhang, et al., 2005). As seen in Fig. 4b, at room temperature, the Hall conductivity σxy reveals clear plateaus at 2e2/h for both electrons and holes, while the longitudinal conductivity ρxx approaches zero. (Quantum Hall effect is measured by σ=υe2h, where “e” is the elementary charge, h is the Planck’ constant, and υ is the “filling factor.” If υ is an integer, it will be an “integer quantum hall effect,” and if υ is a fraction, it will be a “fractional quantum Hall effect.” Here, at room temperature, the filling factor of graphene is υ=2.)


Figure 4.

(a) The resistivity of a single layer of graphene vs. gate voltage. (b) The quantum Hall effect in single-layer graphene. Figures taken from (Novoselov, et al. (2005).

For sensing or transistor application, we should utilize the strong gate dependence of graphene. To do this, we should cut graphene into narrow ribbons because graphene does not have a band gap, and thus resistivity variation is small. However, graphene in narrow ribbons, makes an opening in the band gap that is proportional to the width of the ribbon by increasing the momentum of charge carriers in the traverse direction when shrinking them. This band gap in carbon nanotubes is proportional to its diameter. The opening of a band gap in graphene ribbons has recently been observed in wide ribbon devices lithographically patterned from large graphene flakes (Han, et al., 2007) and in narrow chemically synthesized graphene ribbons (Li, et al., 2008).

3. Molecular Dynamics (MD) modeling

Molecular dynamics is nothing but classical dynamics. Indeed classical dynamic equations of motion are valid for slow and heavy particles, with typical velocities υ<<c, (where c is the speed of light) and masses m>>me, (where me is the electron mass). Therefore, we can use them for atoms, ions, and molecular mass only in slow motion (slower than thermal vibration).

This technique is based on computing potential energy that is typically considered only as a function of the system spatial configuration and is described by means of interatomic potentials. These potentials are considered as known input information; they are either found experimentally or are computed by averaging over the motion of the valence electrons in the ion’s Coulomb field by means of quantum ab initio methods.

The main equation that we utilize in this technique is the Lagrange equations of motion. For a system of N interacting monoatomic molecules, the Lagrange equation turns Newtonian equations divided into “Dissipative equations” and “Generalized Langevin equations.”

Integrals of motion are more functions that are useful for modeling in MD technique. Their notable property is depending only on the initial conditions and staying constant in time. Some of these equations are as follows:

where E is the total energy, P is the total momentum, and M is the total angular momentum in a system with only conservative forces.

4. Lattice Mechanics

Lattice mechanics is an approach to utilize natural symmetries for classical particle mechanics in materials that repetitive or regular atomic structure such as polymers or crystalline materials.

This approach is based on two principles:

Principle 1: Mathematical models such as functions, matrices, operators, and so on are invariant with respect to the symmetry of lattice.

Principle 2: symmetry assumption causes loading is symmetry, and thus we will have symmetry effect.

In this approach, we use the concept of “Unit Cell” instead of particles itself. Unit Cell is an arbitrary part of a lattice, in which we could gain the whole lattice by repeating that in a fixed direction and certain distance.

Hence, we can define displacement vector in a lattice in terms of unit cell as follows:

where n is the primary position of the unit cell:

rn(t), is the position of the unit cell at time t, which is composed of the position of all particles in the unit cell at time t; un(t), is the displacement vector for the unit cell. By calculating the total lattice kinetic energy, we are able to get lattice Lagrangian and consequently the equation of motion as follows:


where U is the total Lattice Potential energy and fnext is the external load on the Unit Cell n. We absolutely insist that this equation is for each unit cell, that is, with changing n we will have a system of motion differential equations.

Using Taylor series for U results in the following:




and n' is indices for neighboring cell.

Since we need neighboring cell in forming equations, we must define another concept as “associate cell”. Associate cell is the smallest part of a lattice that represents its mechanical properties perfectly. In lattice mechanic we restrict our studies to associate cell which can comprise several unit cells; this is the consequence of principle I discussed above.

To solve the above-mentioned equation of motion, we should use mathematics like the Fourier transform or the Laplace transform.

In terms of fext, we have three types of problems:

  1. fext=0 or free lattice: the solution is gained by superposing of normal modes called standing waves, because the lattice oscillate around its equilibrium position.

  2. fext0 : to solve this type of problems we use Green’s function method that requires use of Unit pulse convolution.

  3. Quasi-static problem or approximation: this is a name for time-independent problems where any noticeable change of the external forcing occurs during a period of time that is much longer than the characteristic time of atomic vibrations. This leads to eliminate the first term in motion equation, so we have the following:


To solve this equation, we can use green’s functions method. One of issue problems in quasi-static approximation is multiscale boundary conditions. We discuss this in a separate clause.

5. Multi scale boundary condition

Each method has its own limitations to use. These limitations are only in time and size scale. The multi scale method tries by dividing the model in different areas in terms of size and time scale (coarse-fine grain) and relating them together generates an assimilated method. In other words, we may be able to solve one part of problem with atomistic modeling and the other part with continuum; the multi scale method uses both of them concurrently and then couples them together.

To couple the methods, we define a region of pseudo-atoms called handshake or pad. The position of pseudo-atoms is determined by the finite element method. The handshake has a duty to absorb the fine grain excitation and transfer the effects of coarse grain surrounding boundary.

If u1 is the displacement of the pseudo-atom, u0 is the fine grain displacement, and ua is the coarse grain displacement, then we can say:


where Θ and Ξ are unknown operators.

6. Multi scale Modeling

The philosophy of arising multiple- scale methods is that, in actuality, nano- materials are always used along with large- scale materials and we are compelled to create a method for modeling them. Atomistic methods such as MD and ab initio are not perfect to model the entire configuration because these methods are limited to time and length scales; thus, they are validated only through fine- scale parts of configuration and not for the other part. There is the same situation for continuum methods because they are validated for large time and length scale, and they are not useable for fine- scale parts of these configurations. Here is the point that the role of multiple- scale methods becomes prominent. Multiple scale methods try to blend the methods that are validated on their own scales of time and length separately.

The base of all these methods is that each scale is modeled by its special method, and their output becomes boundary condition for the other; indeed, there will be exchanges between these methods to model entire of the configuration.

Whatever method we take, it must be free of two issue problems:

  • The wavelengths emitted by the MD region are drastically smaller than that can be absorbed by continuum region. It means we have differential energy for these two areas. If our approach is not capable to put in any experience for this redundant energy, it leads to this result that some of the wave will reflect back into the MD domain. This means that we will have spurious heat generation in the MD region and so a mistaken simulation especially in instances of plasticity.

  • Furthermore, the transition from the MD region to the continuum region was followed with extension in timescale, which affects the validity of our relations. Hence, we must take such a method that solves this issue.

6.1. MAAD

MAAD or “macroscopic, atomistic, ab initio dynamics” is a multiple- scale approach that incorporates tight binding (TB), molecular dynamics (MD), and finite elements (FE) concurrently to model a configuration of nano- and large- scale materials. In this method, the FE Mesh will be done until we approach the same size as much as atomic spacing; from here, the MD method is entered until we arrived in a physical phenomenon like a crack tip. At this point, we will use the TB approach. Thus, we will have two overlapping regions: FE/MD and MD/TB. Not only here but also in all multiple- scale methods, such overlapping areas are termed “handshake” regions. In this method, each handshake region provides a contribution to the total energy of the system. This contribution is done by linear law in each handshake region. Thus, the total energy that will be used to find equation of motion is as follows:


There are two problems in this method. First, having finite elements in the scale of atomic space prolongs the process of solving by increasing time steps. Second, that it seems unphysical that continuum relations can be evolved at the same timescales as the atomistic variables.

6.2. Coarse-Grained Molecular Dynamics (CGMD)

This approach couples FE and MD methods together. In this method to derive the governing equations of motion, we use an approximation of coarse-grained energy that converges to the exact atomic energy like the following:



Uint=3(NNnode)kT Internal energy
Mjku˙j.u˙k Kinetic energy Potential energy
u,u˙ Displacement degree of freedom and the velocity

Uint represents the thermal energy of those degrees of freedom of coarse grained that have not been included in FE considered nodes. Obviously, when the number of nodes approaches the number of atoms, this term vanished and the full atomistic energy is recovered.

Finally the equation of motion is then given to be as follows:



Mij Mass matrix
Gij Stiffness like quantity
ηij Time history of memory function
Fi(t) Random force

6.3. Quasi-Continuum Method

This approach is an adaptive FE method. The link between atomistic and continuum is obtained by the use of the Cauchy Born rule. The Cauchy Born rule assumes that the continuum energy density W can be computed using an atomistic potential.

The first Piola-Kirchoff stress tensor in the Cauchy Born rule is

and the Lagrangian tangent stiffness is:


where F is deformation gradient.

The major restriction as well as implication of the Cauchy Born rule is that the deformation of the lattice underlying a continuum point must be homogeneous. This results from the fact that the underlying atomistic system is forced to deform according to the continuum deformation gradient F.

6.4. Coupled Atomistics And Discrete Dislocation (CADD)

This method sets a quasi-static coupling between molecular statics and discrete dislocation plasticity. One of the best assumptions in this approach is that defects within atomistic region are allowed to pass through the atomistic/continuum border into the continuum which is modeled via discrete dislocation mechanics.

Quantities such as stresses, strains, and displacements can be written in the contribution form:

where A is a typical quantity, A˜ is the contribution from the discrete dislocation, and A^ is a correction term we introduce because of the fact that the discrete dislocation solution is for an infinite medium. The continuum energy can be written as

Ec=12Ωc (u^+u˜):(ε^+ε˜)dVdΩT T0(u^+u˜)dA

where T0 is the traction on the continuum boundary dΩΤ.

The equilibrium condition is

Where u˜ is displacement field. Using this equation, we could find displacement fields. By the use of displacement field we are able to find the forces exerted on the discrete dislocation as


At this point, an iterative procedure involving the discrete dislocation positions, FE positions and atomic positions is solved until all degrees of freedom are at equilibrium.

However, we have some challenge in this method as follows:

  • extending this approach to dynamic problems

  • simulating the passage of defect from the atomistic to continuum in three dimensions

  • annular dislocations that reside in both atomistic and continuum regions at the same time

6.5. Bridging Domain

Bridging Domain is a dynamic multiple scale method that couples MD with continuum. In this method the boundary that overlays the MD region and surrounding continuum region has varying size, termed the bridging domain.

Within overlaying bridging domain, the Hamiltonian is defined as a linear combination of the molecular and continuum Hamiltonians:


where α acts to scale the contribution of each domain to the total Hamiltonian.

To make compatibility between the molecular and the continuum regions, we involve Lagrange multipliers in the overlap region as


Where gI are the Lagrange Multipliers, uiJ are the FE nodal displacements, and diJ are the MD displacements.

The coupled equations of motion will be as follows:


Where the standard equations are augmented by the Lagrange multiplier-based constraint forces FIL and fIL. The bar symbols overlaying the FE and MD mass matrices indicate that they need to be modified within the overlapping region.

7. Bridging- Scale Modeling

The bridging scale method that we want to introduce is an incorporation of MD and continuum. The total displacement as our solution is decomposed to fine and coarse- scale as follows:


where u¯ is a coarse- scale solution that will be state in terms of finite element shape functions and u' is a fine- scale solution. These two solutions are orthogonal, that is, one projection onto another is zero.

If x=xα (we use Greek indices for atoms and Roman indices for coarse- scale nodes) was the initial position for typical atoms of interesting body we can rewrite (23) as:


u¯(Xα) is an average behavior of body since it’s interpolation between atoms so we anticipate it is continues function. We can write:


where, NIα=NI(Xα) is a shape function of node I at initial atomic position Xα and dI is the displacement degree of freedom at node I.

As we said, we want to use MD method in this approach, to do this we have to utilize qα. qα is a displacement solution we derive from MD simulation (qα has the same role as u(Xα) in other words, we are simulating u(Xα) by the use of qα. Undoubtedly, qα will have projection onto u¯ (coarse- scale solution) and u (fine- scale solution), and with respect to their orthogonality, each projection will not represent the other one, so we can easily find the fine- scale solution, by subtracting the projection of qα onto coarse- scale solution, from qα, that is,

or with indices representation:


INIαWI is a coarse- scale projection. In other words, to find coarse- scale projection, we have projected qα onto shape functions NI. However, there is one burning question: how can we find WI? The answer is laid in least squares weighted residual method. First, we define J function as follows:


where mα is the atomic mass. The choice of J (called the mass-weighted square of the fine scale) allows the coarse and fine- scale kinetic energies to be decoupled. Minimizing of J will result in WI. If we gather all masses in a matrix as follows:

MA=[m1I3×3      m2I3×3      ]

J can be writhed as:


Solving for w to minimize this quantity gives

where the coarse- scale mass matrix M is

Substituting W in Equation 26 with Equation 31, we will have

where the matrix P has been defined as

So we can write the total displacement uα (with respect to Equation 23-25 and 33) as

The last term in the above expression is called the “bridging scale”; it is that part of the molecular dynamics calculation that must be subtracted from the total in order to achieve a complete separation of scales. Our sake about the “scales” is fine and coarse scales, and our sake about “complete separation” is the orthogonality of these scales.

Indeed, P was the operator for the projection onto the coarse scale; with respect to Equation 35, we could represent another more perfect operator for this projection as



Thus, we can rewrite 35 as

7.1. Equations of motion

In this step, we will derive the coupled MD and FE equations of motions. First, we adopt the Lagrangian equation definition we submitted in MD discussion for multiscale condition as “multiscale Lagrangian.”

7.1.1. Multiscale Lagrangian


Substituting u˙ in Equation 40 with Equation 41 we will have

K(u˙)=12u˙TMAu˙=12d˙TNTMANd˙+d˙TNTMAQq ˙+12q˙TQTMAQq ˙=12d˙TMd ˙+12q˙TM q˙

The second term in the right- hand side of abvoe equation is zero because


The fine- scale mass matrix M in equation 42 defined as


This will be proven by the complete expression for Q (Equation 36) and coarse- scale mass matrix M (Equation 30).

Thus, we can write Lagrangian 36 as

L(d,d˙;q,q˙)=12d˙TMd ˙+12q˙TM q˙U(Nd+Qq)+fextTNd+fextTQq

This is the desired equation for us because we have eliminated cross terms such as d and q in the kinetic energy. The bridging scale is responsible for this elimination. Using the bridging scale makes us capable of decomposing total kinetic energy into the sum of the kinetic energy in the coarse scale plus that in the fine scale.

7.1.2. Multiscale equations of motion

We can derive equations of motion by the use of Lagrangian equation:


Substituting the Lagrangian Equation 44 into Equations 45 and 46 gives us

M q¨=U(d,q)q

For simplicity, we have ignored the external force. Using the chain rule, we can rewrite Equations 47 and 48 as

Md ¨=Uuud
M q¨=Uuuq

As we see, the derivation of potential energy U couples scales in both equations. We know that the variation of energy to displacement is the same type as force. These means we can conclude that the derivation of potential energy is the total forces on atoms:


Using Equations 51 and 38, we can rewrite Equations 49 and 50 as

Md ¨=(ud)TUu=NTf
M q¨=(uq)TUu=QTf

where M =QTMA. If the external forces in Equation 44 are kept in formulation, we will have

Md ¨=NT(f+fext)
M q¨=QT(f+fext)

Equation 55 can be rewritten as


Q is a singular since multiplying a field by Q gives the fine- scale part of the field, and many different fields can have same fine- scale part. Therefore the solution of and equation like Qu=u' for u is nonunique, and Q must be singular. It is concluded that Equation 56 does not have a unique solution, and we are free to choose any q that has the proper fine- scale part and thus satisfies Equation 56. One q that comply these conditions is that which exactly satisfies the following equation:

This is nothing but a molecular dynamics simulation for the atomic displacements, with the atomic forces given by Equation 51. So we are able to solve this equation using MD simulation.

Anyway, now we have the equation of motions as

Md ¨=NT(f+fext)

There are some relevant comments:

  • Equation 58 is a fine- scale equation of motion and a standard MD solver can be used to obtain the MD displacements q.

  • Equation 59 is a coarse- scale equation of motion. This is an FE momentum equation, and to solve it, we should use finite elements method. It should be noted that for the coarse scale, the summation sign turns into integral sign.

  • As we said before, q (MD solution) simulate u (real total displacement). Although we proved Equation 58 for q, it is obvious that u as real total displacement have to satisfy Equation 58 in each atomic position (Xα), this is found from Newton’s second law. u is not necessarily equal to q; they just have the same behavior like solutions of a differential equations. Indeed, initial conditions determine their equality. Thus, if these two quantities have the same initial conditions, then (which they should) then they are equal forever.

  • Equality of q and u and 38 gives

It means that the coarse scale is a projection of q, and instead of solving the FE equation, we can use q. However, as we will discuss, since we will eliminate the fine scale from the coarse scale, q will be limited to fine scale (not the entire domain), and thus it is not possible to calculate the coarse- scale solution everywhere via direct the projection of the MD displacements.

  • Equations 58 and 59 converge to each other if FE nodal positions approach MD atomic positions.

  • This convergence makes us to conclude FE equation (Equation 59) overlays both fine and coarse scales. This is the unwanted degree of freedom that prolongs spent time for FE solution, so we should eliminate the effects of the fine scales that lie in coarse- scale region.

7.2. Removing fine- scale degree of freedom in coarse- scale region

In the previous section, we saw the fine- scale equation of motion (Equation 58) and coarse- scale equation of motion (Equation 59) in areas among this two- scale overlay together, and this is in contrast with the orthogonality of these two scale, so we should eliminate the effect of fine scale from coarse scale. However why do we eliminate fine scale from coarse scale and not vice versa? The answer is because eliminating one scale from the other will enter another terms in the eliminated scale’s equation of motion. Since we want to avoid wave reflection back into the MD domain and avoid heat retention within this region, we should add some terms to the fine- scale equation of motion, and it is the best opportunity for us to achieve this gold by eliminating the fine- scale effect from the coarse- scale domain. Thus, we limit our studies to the MD region and its equation of motion (Equation 58) with assumption fext=0:

To find the area that shall be removed, we should consider the behavior of these two scales. In fine scale, we have a harmonic or nonlinear potential in atomic interactions. However, at coarse scale, we have a harmonic or linear potential energy in atomic interaction. Thus, we can conclude that there is a transition area among these two scales that belong in the MD region. This is the area that shall be eliminated, and instead of this elimination, we should add some boundary conditions for the MD region. Thus, we divide the MD region in two sections: linear and nonlinear. This division is called the “linearized MD equation of motion.” Two find MD boundary condition, since we are transient from fine to coarse scale in the transition area, we use FE method and then Lattice mechanics. This means we are treating with transition area, like a coarse scale, and maybe we could say we are dividing the MD area into fine and coarse sections. We first study a chain of atoms and then consider the general state.

Consider a one-dimensional chain of atoms according to the following figure:


From the lattice mechanics, we have


where fn denotes the total force on atom (unit cell) n and Knn' is the springy factor. Using the Lagrange equation results in the equation of motion as follows:


This equation is validated only for harmonic lattice, and fnext(t) is the total external force acting on the unit cell n.

By dividing the MD region into fine and coarse scales, we have

where u¯ denotes the coarse scale and u denotes the fine scale of the MD region.

For the left side of Equation 61, we can write




Moreover, fext is the total external force acting on the harmonic section of the MD by an un-harmonic section of the MD (note that fext is different from fext we neglected in Equation 48; fext takes place within the MD section while fext lies outside of this section).

A comparison between Equations 66 and 61 gives


In Equation 66, we have grouped terms based on fine scale and coarse scale. Because of the orthogonality of coarse and fine scales and that the timescale for the coarse scale is much larger than that of atomic vibrations present in the fine scale, it is reasonable to equate the corresponding parts from each side of equality sign, that is,


Equation 70 is the one we were looking for; in other words, by removing the fine- scale degree Of freedom in the coarse- scale region, Equation 61 turns into Equation 70. Thus we use Equation 70 instead of Equation 61. However we still have not terminated; we should calculate fext as an issue problem. To do this, we use lattice mechanics.

In lattice mechanics, Equation 70 turns into Equation 71:


Equation 70 is the one we were looking for; in other words, by removing the fine- scale degree Of freedom in the coarse- scale region, Equation 61 turns into Equation 70. Thus we use Equation 70 instead of Equation 61. However we still have not terminated; we should calculate fext as an issue problem. To do this, we use lattice mechanics.

In lattice mechanics, Equation 70 turns into Equation 71:


where un is the displacement of unit cell n.

Again, consider the one- dimensional chain of atoms like the following figure:


n=0 is the position that elimination has took place, and the result of that is fnext. We can write the following:


Applying discrete Fourier transform (DFT) to Equation 71, we have


Using Laplace transform (LT) for Equation 73 and solving for resulting displacements U^ gives






For f0ext we can write:


Taking LT from Equation 77 gives impedance force:


So to find f0ext, we should calculate U1´(s), This requires to omit F0ext(s) from Equation 74, so we first take the inverse discrete Fourier transform of Equation 61 in terms of p:




and N is the total number of unit cells in the domain. Now we write Equation 66 for n=0,1:


Solving these two equations simultaneously can omit F0ext(s) and gives


We can rewrite Equation 83 as




Now taking the inverse Laplace transform (ILT) of Equations 78 and 84 and convoluting the displacement, we obtain


We know


And then


Where the time history kernel θ(t) and random force R0f(t) are defined to be


Other terms in Equation 71 is equal with f(u), and we can rewrite Equation 71 as


We can turn u to q as


Thus, the final coupled forms of the equations of motion are


where the impedance force f0imp acts only on the boundary atom 0.

7.2.1. Elimination of fine- scale degrees of freedom: 3D generalization

In 3D state, each unit cell can be labeled with three indices, l, m, and n indicating the position along axes in the direction of the three primitive vectors of the crystal structure. The equation of motion can be rewritten as


where μ and υ represent the range of the forces in the m and n coordinate directions. We consider the one dimensional boundary condition in our study. Our boundary is located in l=0, like the following figure:


Like the one-dimensional chain, we can write


By taking LT and DFT of Equation 96 and solving for U^', we obtain






For f0,m,next, we can write


Taking LT from Equation 101 gives impedance force:


Thus, to find f0,m,next we should calculate U'1,m',n'(s). This requires to omit F^0ext from Equation 74, so we first take inverse discrete Fourier transform of (98) in terms of p:

U˜l'(q,r,s)=MA1G ˜l(q,r,s)F^0ext(q,r,s)+R ˜ld(q,r,s)


R ˜ld(q,r,s)=l'=L/2L21G ˜ll'(q,r,s)(sUl''(q,r,0)+U˙l''(q,r,0))

Here L is the total number of unit cells in the x-direction. Letting l=0,1 and eliminating F^0ext we have:

U˜'1(q,r,s)=Q˜(q,r,s)(U˜0'(q,r,s)R ˜0d(q,r,s))+R ˜1d(q,r,s)


Q ˜(q,r,s)=G ˜1(q,r,s)G ˜01(q,r,s)

Taking IFT of 105 and using convolution property of the DFT, we gain


where upon plane l we can write:


Taking ILT of Equation 108 gives

Rl,m,nd(t)=l'=L/2L21m'=M/2M21n'=N/2N21(g ˙ll',mm',nn'(t)ul',m',n''(0)+gll',mm',nn'(t)u˙l',m',n''(0))


g ˙l,m,n(t)=L1(sGl,m,n(s))

Substituting Equation 107 into Equation 102 and applying ILT, the impedance force boundary condition obtains as


where the time history kernel θ(t) is defined to be




where Rd is given by Equation 109.

For simplicity, we rewrite Equation 112 as


where nc refers to a maximum number of atomic neighbors, which will be used to compute the impedance force. Certainly, all of this nc atoms are located in the MD region with linear or harmonic behavior.

As we evolved in Equation 92, we can write the final coupled form of the equations of motion as:




Knowing that



  1. The impedance force fimp became an appliance for us to vanish all fine- scale information that cannot be represented in the continuum. fimp=0 means we do not have any unharmonic area in the MD region.

  2. Rf and Rd are called random or stochastic forces. These forces arise from the initial conditions in the eliminated fine- scale degrees of freedom. These terms represent thermally motivated displacements and forces exerted on the reduced MD system by the removed fine- scale degrees of freedom.


1 - Allen, Matthew J., Vincent C. Tung, and Richard B. Kaner. "Honeycomb carbon: a review of graphene." Chemical reviews 110.1 (2009): 132 145.
2 - Geim, Andre Konstantin. "Graphene: status and prospects." science 324.5934 (2009): 1530 1534.
3 - Rao, C. N. R, et al. "Graphene: The New Two‐Dimensional Nanomaterial." Angewandte Chemie International Edition 48.42 (2009): 7752-7777.
4 - Stankovich, Sasha, et al. "Synthesis and exfoliation of isocyanate-treated graphene oxide nanoplatelets." Carbon 44.15 (2006): 3342 3347.
5 - Stankovich, Sasha, et al. "Synthesis of graphene-based nanosheets via chemical reduction of exfoliated graphite oxide." Carbon 45.7 (2007): 1558 1565.
6 - Ferrari, Andrea C. "Raman spectroscopy of graphene and graphite: disorder, electron–phonon coupling, doping and nonadiabatic effects." Solid state communications 143.1 (2007): 47 57.
7 - Katsnelson, M. I., and K. S. Novoselov. "Graphene: new Bridge between condensed matter physics and quantum electrodynamics." Solid State Communications 143.1 (2007): 3 13.
8 - Bunch, J. Scott, et al. "Electromechanical resonators from graphene sheets." Science 315.5811 (2007): 490 493.
9 - Atalaya, Juan, Andreas Isacsson, and Jari M. Kinaret. "Continuum elastic modeling of graphene resonators." Nano Letters 8.12 (2008): 4196 4200.
10 - Sakhaee-Pour, A., M. T. Ahmadian, and A. Vafai. "Applications of single-layered graphene sheets as mass sensors and atomistic dust detectors." Solid State Communications 145.4 (2008): 168 172.
11 - Ramanathan, T., et al. "Functionalized graphene sheets for polymer nanocomposites." Nature nanotechnology 3.6 (2008): 327 331.
12 - Saremi, F., et al. "Adsorption of carbon monoxide on a (6, 6) armchair carbon nanotube: ab initio study." Journal of Physical and Theoretical Chemistry Islamic Azad University of Iran 4.4 (2008): 235 238.
13 - Ekinci, K. L. "Electromechanical transducers at the nanoscale: actuation and sensing of motion in nanoelectromechanical systems (NEMS)." small 1.8 9 (2005): 786 797.
14 - Eom, Kilho, et al. "Nanomechanical resonators and their applications in biological/chemical detection: nanomechanics principles." Physics Reports 503.4 (2011): 115 163.
15 - Arlett, J. L., E. B. Myers, and M. L. Roukes. "Comparative advantages of mechanical biosensors." Nature nanotechnology 6.4 (2011): 203 215.
16 - Poot, Menno, and Herre SJ van der Zant. "Mechanical systems in the quantum regime." Physics Reports 511.5 (2012): 273 335.
17 - Peng, H. B., et al. "Ultrahigh frequency nanotube resonators." Physical Review Letters 97.8 (2006): 087203.
18 - Sakhaee-Pour, A., M. T. Ahmadian, and A. Vafai. "Potential application of single-layered graphene sheet as strain sensor." Solid State Communications 147.7 (2008): 336 340.
19 - Sheehan, Madoc E., and Paul N. Sharratt. "Molecular dynamics methodology for the study of the solvent effects on a concentrated Diels-Alder reaction and the separation of the post-reaction mixture." Computers and chemical engineering 22 (1998): S27-S33.
20 - Kresse, Georg, and Jürgen Hafner. "Ab initio molecular dynamics for liquid metals." Physical Review B 47.1 (1993): 558.
21 - Gumbsch, Peter, and Rowland M. Cannon. "Atomistic aspects of brittle fracture." MRS Bulletin 25.05 (2000): 15 20.
22 - Müser, Martin H. "Towards an atomistic understanding of solid friction by computer simulations." Computer physics communications 146.1 (2002): 54 62.
23 - Horstemeyer, M. F., M. I. Baskes, and S. J. Plimpton. "Computational nanoscale plasticity simulations using embedded atom potentials." Theoretical and applied fracture mechanics 37.1 (2001): 49 98.
24 - Abraham, Farid F., et al. "Spanning the continuum to quantum length scales in a dynamic simulation of brittle fracture." Europhysics Letters 44.6 (1998): 783.
25 - Broughton, Jeremy Q., et al. "Concurrent coupling of length scales: methodology and application." Physical review B 60.4 (1999): 2391.
26 - Tadmor, Ellad B., Michael Ortiz, and Rob Phillips. "Quasicontinuum analysis of defects in solids." Philosophical magazine A 73.6 (1996): 1529 1563.
27 - Tadmor, E. B., Rob Phillips, and M. Ortiz. "Mixed atomistic and continuum models of deformation in solids." Langmuir 12.19 (1996): 4529 4534.
28 - Shenoy, V. B., et al. "Quasicontinuum models of interfacial structure and deformation." Physical Review Letters 80.4 (1998): 742.
29 - Miller, R., et al. "Quasicontinuum models of fracture and plasticity." Engineering Fracture Mechanics 61.3 (1998): 427 444.
30 - Miller, R., et al. "Quasicontinuum simulation of fracture at the atomic scale." Modelling and Simulation in Materials Science and Engineering 6.5 (1998): 607.
31 - Ortiz, M., et al. "Mixed Atomistic–Continuum models of material behavior: The art of transcending atomistics and informing continua." Mrs Bulletin 26.03 (2001): 216 221.
32 - Wagner, Gregory J., and Wing Kam Liu. "Coupling of atomistic and continuum simulations using a bridging scale decomposition." Journal of Computational Physics 190.1 (2003): 249 274.
33 - Geim, Andre Konstantin. "Graphene: status and prospects." science 324.5934 (2009): 1530 1534.
34 - Hirth, John P., and Jens Lothe. "Theory of dislocations." (1982).
35 - Sholl, David S. "Structural Nanocrystalline Materials. Fundamentals and Applications. Von Carl C. Koch, Ilya A. Ovid'ko, Sudipta Seal und Stan Veprek." Angewandte Chemie 120.6 (2008): 1022 1022.
36 - Zhu, T., and Ju Li. "J. Ultra-strength materials Prog." Mater. Sci 55 (2010): 710 757.
37 - Lee, Changgu, et al. "Measurement of the elastic properties and intrinsic strength of monolayer graphene." science 321.5887 (2008): 385 388.
38 - Liu, Fang, Pingbing Ming, and Ju Li. "Ab initio calculation of ideal strength and phonon instability of graphene under tension." Physical Review B 76.6 (2007): 064120.
39 - Kim, Kwanpyo, et al. "Ripping graphene: preferred directions." Nano letters 12.1 (2011): 293 297.
40 - Sen, Dipanjan, et al. "Tearing graphene sheets from adhesive substrates produces tapered nanoribbons." Small 6.10 (2010): 1108 1116.
41 - Meyer, Jannik C., et al. "Direct imaging of lattice atoms and topological defects in graphene membranes." Nano letters 8.11 (2008): 3582 3586.
42 - Warner, Jamie H., et al. "Dislocation-driven deformations in graphene." Science 337.6091 (2012): 209 212.
43 - Kim, Kwanpyo, et al. "Grain boundary mapping in polycrystalline graphene." ACS nano 5.3 (2011): 2142 2146.
44 - Huang, P. Y., et al. "Grains and grain boundaries in single-layer graphene atomic patchwork quilts." Nature 469 (2011): 389 392.
45 - Ruiz-Vargas, Carlos S., et al. "Softened elastic response and unzipping in chemical vapor deposition graphene membranes." Nano Letters 11.6 (2011): 2259 2263.
46 - Warner, Jamie H., et al. "Dislocation-driven deformations in graphene." Science 337.6091 (2012): 209 212.
47 - Ong, Zhun-Yong, Eric Pop, and Junichiro Shiomi. "Reduction of phonon lifetimes and thermal conductivity of a carbon nanotube on amorphous silica." Physical Review B 84.16 (2011): 165418.
48 - Kang, Kwangu, et al. "Lifetimes of optical phonons in graphene and graphite by time-resolved incoherent anti-Stokes Raman scattering." Physical Review B 81.16 (2010): 165405.
49 - Qiu, Bo, and Xiulin Ruan. "Reduction of spectral phonon relaxation times from suspended to supported graphene." Applied Physics Letters 100.19 (2012): 193101.
50 - Tohei, Tetsuya, et al. "Debye temperature and stiffness of carbon and boron nitride polymorphs from first principles calculations." Physical Review B 73.6 (2006): 064304.
51 - Nicklow, R., N. Wakabayashi, and H. G. Smith. "Lattice dynamics of pyrolytic graphite." Physical Review B 5.12 (1972): 4951.
52 - Nihira, Takeshi, and Tadao Iwata. "Temperature dependence of lattice vibrations and analysis of the specific heat of graphite." Physical Review B 68.13 (2003): 134305.
53 - Chen, Shanshan, et al. "Raman measurements of thermal transport in suspended monolayer graphene of variable sizes in vacuum and gaseous environments." ACS nano 5.1 (2010): 321 328.
54 - Balandin, Alexander A. "Thermal properties of graphene and nanostructured carbon materials." Nature materials 10.8 (2011): 569 581.
55 - Chen, Shanshan, et al. "Thermal conductivity of isotopically modified graphene." Nature materials 11.3 (2012): 203 207.
56 - Pettes, Michael Thompson, et al. "Influence of polymeric residue on the thermal conductivity of suspended bilayer graphene." Nano Letters 11.3 (2011): 1195 1200.
57 - Anthony, T. R., et al. "Thermal diffusivity of isotopically enriched C 12 diamond." Physical Review B 42.2 (1990): 1104.
58 - Berman, R. "Thermal conductivity of isotopically enriched diamonds." Physical Review B 45.10 (1992): 5726.
59 - Pop, Eric. "Energy dissipation and transport in nanoscale devices." Nano Research 3.3 (2010): 147 169.
60 - Chen, Z., et al. "Thermal contact resistance between graphene and silicon dioxide." Applied Physics Letters 95.16 (2009): 161910.
61 - Koh, Yee Kan, et al. "Heat conduction across monolayer and few-layer graphenes." Nano Letters 10.11 (2010): 4363 4368.
62 - Mak, Kin Fai, Chun Hung Lui, and Tony F. Heinz. "Measurement of the thermal conductance of the graphene/SiO2 interface." Applied Physics Letters 97.22 (2010): 221904 221904.
63 - Liao, Albert D., et al. "Thermally limited current carrying ability of graphene nanoribbons." Physical Review Letters 106.25 (2011): 256801.