Strain gauge fields for rippled graphene membranes under central mechanical load: an approach beyond first-order continuum elasticity
We study the electronic properties of rippled freestanding graphene membranes under central load from a sharp tip. To that end, we develop a gauge field theory on a honeycomb lattice valid beyond the continuum theory. Based on the proper phase conjugation of the tight-binding pseudospin Hamiltonian, we develop a method to determine conditions under which continuum elasticity can be used to extract gauge fields from strain. Along the way, we resolve a recent controversy on the theory of strain engineering in graphene: There are no K-point dependent gauge fields. We combine this lattice gauge field theory with atomistic calculations and find that for moderate load, the rippled graphene membranes conform to the extruding tip without significant increase of elastic energy. Mechanical strain is created on a membrane only after a certain amount of load is exerted. In addition, we find that the deformation potential –even when partially screened– induces qualitative changes on the electronic spectra, with Landau levels giving way to equally-spaced peaks.
The interplay of electronic and mechanical properties of graphene membranes is a subject under intense experimental and theoretical investigation Meyer et al. (2007); Bunch et al. (2007); Pereira and Castro-Neto (2009); Guinea et al. (2010); Kitt et al. (2012); Suzuura and Ando (2002); Vozmediano et al. (2010); (8); de Juan et al. (2012a). Mechanical strain induces gauge fields in graphene that affect the dynamics of charge carriers Suzuura and Ando (2002); Pereira and Castro-Neto (2009); Guinea et al. (2010); Vozmediano et al. (2010); Kitt et al. (2012). As graphene can sustain elastic deformations as large as 20% Kim et al. (2009), the resulting pseudo-magnetic fields are much larger than those magnetic fields available in state-of-the-art experimental facilities (for example, the highest magnetic field created at the US National High Magnetic Field Laboratory is slightly larger than 100 Tesla). The presence of a pseudo-magnetic field is observed via broad Landau levels (LLs) in strained graphene nanobubbles on a metal substrate Levy et al. (2010). In addition to the pseudo-magnetic vector potential , strain also induces a scalar deformation potential Suzuura and Ando (2002); de Juan et al. (2007); Choi et al. (2010) that affects the electron dynamics in complex ways.
The theoretical formalism has been laid out within the context of first-order continuum elasticity Suzuura and Ando (2002); Pereira and Castro-Neto (2009); Guinea et al. (2010); Vozmediano et al. (2010); Kitt et al. (2012); de Juan et al. (2007); (8); de Juan et al. (2012a). It is possible to improve the theory from a mechanical perspective. Such a development, provided on the present manuscript, improves our physical understanding of the inter-relation between mechanics and the electrons in graphene. The purpose of the present paper is twofold: First we motivate, build, and validate a novel framework to lay out a theory valid beyond first-order continuum mechanics. This novel formulation brings to the spotlight some of the inherent assumptions of the prevailing theoretical framework; assumptions that have remained to some extent hidden within the continuum formalism. We disclose upfront that the formalism does not take into account the effects of curvature within the framework of Refs. Vozmediano et al. (2010); de Juan et al. (2007); (8); de Juan et al. (2012a); we will address such shortcoming in the near future. Nevertheless, the reader will realize that the inherent formulation of the theory on the present paper remains novel, bringing a deeper understanding of the formalism for studying the effects of mechanical strain on the electronic properties of graphene.
Following recent experimental developments in which graphene membranes are studied with local scanning tunneling microscopy probes Xu et al. (2012); Zan et al. (2012); Klimov et al. (2012), our second goal is to demonstrate this novel formalism on freestanding graphene membranes under load by a sharp tip. The input for this formalism is direct ‘raw’ atomic displacements upon strain, as opposed to the always present continuum deformation field Meyer et al. (2007); Bunch et al. (2007); Pereira and Castro-Neto (2009); Guinea et al. (2010); Kitt et al. (2012); Vozmediano et al. (2010); de Juan et al. (2012a).
The presentation is given in modular and self-contained form. Hence, discussion of the formalism is given first, then the pure mechanics of freestanding membranes is presented, and predictions from the formalism as pertains to freestanding membranes follows. This helps in focusing either on the basic formulation, or on the predictions from this theory on a experimentally-relevant system. Conclusions are given at the end of the manuscript.
Ii What are the underlying assumptions of the theory?
In order to motivate the developments presented here, we express in an explicit form the underlying assumptions of the theory, which can be found as opening statements in Ref. Guinea et al., 2010: “If a mechanical strain varies smoothly on the scale of interatomic distances, it does not break sublattice symmetry but rather deforms the Brillouin zone in such a way that the Dirac cones located in graphene at points and are shifted in opposite directions.” (See also Ref. Castro-Neto et al., 2009.)
Previous statement tells us that –provided strain preserves sublattice symmetry– one can understand the effects of mechanical strain on the electronic structure in terms of a semiclassical approach, as follows: The local strain-induced fields and are incorporated into a spatially-varying pseudospin Hamiltonian , where is the low-energy expansion of the Hamiltonian in reciprocal space in the absence of strain. We will mention a number of times that the semiclassical approximation is justified if the strain is slowly varying, that is, when it extends over many unit cells Suzuura and Ando (2002) and preserves sublattice symmetry Guinea et al. (2010); Vozmediano et al. (2010).
Evidently, it is also possible to determine the electronic properties directly from a tight-binding Hamiltonian in real space, without resorting to the semiclassical approximation and without imposing an a priori lattice symmetry. That is, while the semiclassical is defined in reciprocal space (thus assuming some reasonable preservation of crystalline order), the tight-binding Hamiltonian in real space is more general and can be used for membranes with arbitrary spatial distribution and magnitude of the strain.
Constituting one of the main arguments of the present paper, we show how to determine if mechanical distortions preserve the fundamental sublattice symmetry. We do this by computing, at each unit cell, angular and length changes from the adequate nearest-neighbor vectors. Such measures will become relevant for the strongly inhomogeneous strain created by local probes Xu et al. (2012); Zan et al. (2012); Klimov et al. (2012), and will give a quantitative meaning –for the first time– to statements such as “long-range mechanical distortionSuzuura and Ando (2002)” and “mechanical distortions preserving sublattice symmetry Guinea et al. (2010).” The program comes down to re-expressing the theory beyond continuum elasticity and explicitly on the atomic lattice, such that matters of spatial scale can be analyzed. For clarity we say that –as a matter of definition– there is no explicit information of interatomic distances on a continuum media, so sublattice symmetry cannot be determined on this formulation of the theory.
Indeed, in the only known formulation of the theory (commonly referred to as the lattice, or tight-binding approach), both and are expressed in terms of a continuous displacement field obtained within first-order continuum elasticity (CE) Suzuura and Ando (2002); Pereira and Castro-Neto (2009); Castro-Neto et al. (2009); Guinea et al. (2010); Vozmediano et al. (2010). It is not possible to assess sublattice symmetry on a continuum media, and therefore proper phase conjugation of pseudospin Hamiltonians becomes an implicit assumption of the theory. Continuum elasticity is based on the fundamental assumption, known as Cauchy-Born rule (CBR), that deformations around any material point are homogeneous. But CBR does not hold exactly on the honeycomb lattice, nor under central load or rippling (15).
As an additional contribution on the present paper that adds physical value to our formulation, we mention that an argument was made in the recent past for the inclusion of additional K-point dependent terms to pseudo-magnetic fields Kitt et al. (2012). Working directly on the atomic lattice, it is easy to show that such terms vanish to first order. We will also show how the formalism based on CE Suzuura and Ando (2002); Pereira and Castro-Neto (2009); Castro-Neto et al. (2009); Guinea et al. (2010); Vozmediano et al. (2010) becomes a limiting case of the one presented here, when the distortion at all unit cells is small in comparison to the lattice constant .
Iii Formulating a theory beyond first-order continuum elasticity
iii.1 Sublattice symmetry and measures for long-range mechanical strain
Consider the tight-binding Hamiltonian in reciprocal space with no strain Castro-Neto et al. (2009):
with eV. The relevant vectors are shown in Fig. 1. Note that in this Figure the zigzag direction lies along the y-axis (a more common choice Guinea et al. (2010); Vozmediano et al. (2010) is to have the zigzag direction parallel to the x-axis; this is a minor detail, that has to be kept in mind when comparing our final expressions for gauge fields to previous ones Guinea et al. (2010); Vozmediano et al. (2010).) With the choices for the lattice vectors made in Fig. 1(a) we have , . To test the purported K-point dependent correction to the theory Kitt et al. (2012), we write down all six K-points explicitly:
It follows that:
and . The low-energy expansion of Eqn. (1) (when ) expresses the dynamics of pseudospinor on the honeycomb lattice; in that limit . Though certainly redundant due to crystal symmetry in the absence of strain, one is free to define one at each unit cell, with the finite number () of pseudospin Hamiltonians for a membrane with a finite number () of atoms.
This finite number of pseudospinor Hamiltonians that can be defined on a membrane with atoms represents the first departure of our formulation of the theory when compared with the formalism developed on a continuum media: In the continuum approach Suzuura and Ando (2002); Pereira and Castro-Neto (2009); Guinea et al. (2010); Vozmediano et al. (2010); Kitt et al. (2012); de Juan et al. (2007); (8); de Juan et al. (2012a), the two degrees of freedom of a given local pseudospin Hamiltonian should correspond to those of an underlying unit cell with two atoms. However, the pseudospin Hamiltonian is defined as a continuous function of coordinates and is hence detached from the actual spatial structure of the lattice. This semiclassical approach is justified when the spatial variation of the strain is small on the scale of the lattice constant . In the present work, we develop a more general method which preserves the spatial scale of mechanical distortion relative to , as well as the total number of local pseudospin Hamiltonians. By doing so we can analyze situations in which the continuum approach breaks down. In addition, we show that pseudo-magnetic vector fields do not depend on points.
The only way to know whether the strain preserves sublattice symmetry Guinea et al. (2010) is by analyzing relative atomic displacements. The nearest neighbor vectors for atom () become: , , and (, , and ); see Fig. 1(b). While by construction, is not necessarily equal to for arbitrary strain. To better quantify the local departures from the sublattice symmetry at any given unit cell, we define the differences in angular orientation and length for nearest-neighbor vectors under mechanical load: Writing and for :
where is a unit vector along the z-axis, and:
We reiterate that the existing theorySuzuura and Ando (2002); Guinea et al. (2010); Vozmediano et al. (2010) requires sublattice symmetry to hold: , and . In practice however, as no measure existed to test those requirements, in applying the theory one actually is lead to assume a priori that , and .
Later on, we will have the opportunity to see how quantifiable deviations of sublattice symmetry occur under central load. Spatial locations where and are much larger than zero indicate that the continuum theory ceases to be applicable there, as the lack of sublattice symmetry will not allow proper phase conjugation of pseudospin Hamiltonians. This should not be even surprising because for a reciprocal space to exist one has to preserve the crystal symmetry. When the crystal symmetry is strongly perturbed, the reciprocal space representation looses its physical meaning. In such scenario (and hence the LDOS) is still meaningful, and so is , but starts to be ill-defined as a local lack of sublattice symmetry necessarily implies the lack of proper phase conjugation. Unfortunately, it is not always possible to express distortions using a continuum theory for lattices with inner structure (sublattices and ): first-order CE breaks down (15), and inhomogeneity of the atomic displacements –reflecting lack of periodicity upon non-uniform strain and shown in Fig. 1– sets in.
iii.2 Relative shift of the and points upon strain
In the more general and lattice-explicit approach being presented here, can be obtained at unit cells in which and , by a (local) replacement of in Eqn. (1) with displaced vectors at each of the two sublattice atoms. One realizes that under load the lattice vectors become ; , which in turn leads to renormalized points. So, to first order in displacements, the reciprocal lattice vectors can be obtained from:
where , and . The form of is convenient, as then and therefore using Eqn. (2) one gets:
so that “the Dirac cones located in graphene at points and are shifted in opposite directions Guinea et al. (2010); Castro-Neto et al. (2009).” This fact builds onto the consistency of the present formulation of the theory.
iii.3 Lattice-explicit strain-gauge potentials (negligible curvature)
How does the pseudospinor Hamiltonians look in the new formalism? Standard manipulation ( ; ; , ) leads to the off-diagonal term:
Here, is the change of the hopping parameter upon strain. The other off-diagonal term is:
(note that in Eqn. 11, is expressed in terms of unprimed ’s). When it follows that . If the sublattice symmetry does not hold to measurable extent, Eqn. (III.3) would not be exactly conjugated to Eqn (11), and applicability of the theory at those unit cells is questionable.
Remarkably, the term:
leads to the linear dispersion because the phasors on add up to zero (this can be shown by explicit calculation). Neglect of the term linear on in Ref. Kitt et al. (2012) led to artificial gauges. Hence, it is just the single term:
that leads to pseudo-magnetic gauge field to lowest-order.
When the zigzag direction is parallel to the axis Guinea et al. (2010); Vozmediano et al. (2010), the real (imaginary) part of Eqn. 14 directly leads to the x- (y-)component of . With the choice made in Fig. 1(a) one obtains reverted components:
with the flux quantum. The ‘+’ sign appears for , and ; the ’’ sign (implying ) near , and . The net pseudo-magnetic field is zero. Dirac’s eqn. in terms of is ():
with the identity, and (Eqn. 3.7 in Ref. Suzuura and Ando (2002), or Eqn. 56 in Ref. Vozmediano et al. (2010)). Suzuura and Ando (2002); Guinea et al. (2010); Vozmediano et al. (2010), hence we arrive at:
so at a given unit cell, each component of takes a single value. We assume to be linearly-dependent to the average bond increase Choi et al. (2010):
[ in Fig. 7(a) is similar to the profile in Ref. de Juan et al. (2007).] Eqns. (17) and (18) express the gauge fields in terms of lattice displacements, representing one of our main results. Besides the inherent physical motivation which has been explained in detail, Eqns. (17) and (18) obviate the need for a continuous deformation field, and hold regardless of the magnitude of the deformation, even in the anharmonic regime (refer to Fig. 4(a)).
iii.4 Limiting form of the vector potential for strain varying slowly with respect to
In the limit the theory from CE is restored. Indeed,
and after simple algebraic manipulations one gets:
The novel formalism has been completely motivated, laid out, and validated at this moment. Now, to plot a flattening procedure and a method of finite differences were developed so that the three-dimensional ’s from atomic displacements could be used in Eqns. (16) and (17). See Figure 2(b).
The idea is to flatten locally the three nearest neighbor vectors from their positions under stress such that their lengths and relative angles are preserved to the greatest extent possible. Flattened vectors are then used in the two-dimensional . It should be clear that this procedure can only be appropriate when the vectors () are small with respect to the lattice constant. Extraction of is only sensible at locations away from the tip, where . The process involves rotating the three vectors first so that the three outer points defining these vectors lie on the x-y plane. This is described as step (i-iii) in Figure 2. Once these points lie on the same plane, we perform an additional rotation about the z-axis so that one of the vectors is close to its original projection along the x-y plane (Figure 2(iv)). The process is completed by bringing the atom in the center towards the x-y plane, by setting its magnitude along the z-axis to be zero (Figure 2(v)). The vertical displacement in Fig. 2(iv) is exaggerated (more below).
As said before, is obtained in terms of finite differences on the flattened membrane, as follows:
with a unit vector pointing out of plane. and are to be computed at unit cells for which sublattice symmetry is reasonably preserved. The partial derivative is estimated after flattening as:
See Figure 2(b) for a schematic illustration of the locations involved. The computation of the curl in terms of finite differences reflects the inherently discrete nature of the present formulation. No other result, including the LDOS, required flattening. We concur that the present method for obtaining could benefit from the ideas within the geometrical approach Vozmediano et al. (2010); (8); de Juan et al. (2012a). We expect to address this aspect in the near future. It should be clear, nevertheless, that our approach beyond continuum elasticity never looses its novelty and value.
Motivated by recent experiments Xu et al. (2012); Zan et al. (2012); Klimov et al. (2012), we illustrate previous considerations by studying the freestanding graphene membranes under central load by a sharp Scanning Tunneling Microscope (STM) tip. As we will show, the membranes are rippled before load because of dynamic (temperature-induced) structural distortions Fasolino et al. (2007), and because of static structural distortions created by interaction with a substrate, the deposition process Meyer et al. (2007), or line stress at edges.
Iv The mechanical behavior of Freestanding graphene membranes under central mechanical load
iv.1 Details of the systems studied and the molecular dynamics calculations
We considered triangular membranes with side nm and 0.16 million atoms [Fig. 3(a)]. We have chosen triangular boundaries since they are known to create the most uniform pseudo-magnetic field Guinea et al. (2010). Equilibrium atomic configurations were obtained from classical molecular dynamics simulations at 1 Kelvin Plimpton (1995). Prior to load, the initially flat membrane is allowed to relieve line strain at its edges, equilibrating forces for 500,000 fs with all atoms moving freely. At this low temperature the average lattice constant is Å. In the initial state after relaxation, the membrane is rippled, with a minimum-to-maximum vertical displacement of 0.8 nm (see leftmost subplot in Fig. 3(c)).
The membrane rims (shown in brown in Fig. 3(a)) represent the mechanical support of a freestanding membrane; they are clamped after the equilibrium rippled conformation is obtained. The height fluctuations seen on the first subplot in Fig. 3(c) tell us that a finite-size graphene membrane behaves as a shell in equilibrium, because it has nonzero local curvature in the absence of applied strain. This behavior will be necessarily linked to the magnitude of the mechanical strain upon load. The (static) rippling discussed here and due to finite size is different from the dynamic effect produced by temperature Fasolino et al. (2007). We must note that most theoretical works consider as their starting point a planar membrane (a thin plate in mechanical jargon). Exceptions are presented in the geometrical approach (Refs. de Juan et al. (2007); Vozmediano et al. (2010); de Juan et al. (2012a)), a formulation of the theory still on a continuum media, where higher-order terms –related to curvature– enter in. Being a theory on a continuum as well, the issue of the scale of the mechanical distortion here discussed carries on.
Strain is induced on the rippled membrane by pushing down a spherical tip (3 nm in diameter), interacting with the membrane via a van der Waals term (details can be provided upon request). The tip pushes the membrane at speed nm/fs to a distance , where is the load time. The load protocol used here is different than the one used in experiments Xu et al. (2012); Klimov et al. (2012), where the tip retracts away from the membrane. The membrane –initially 0.2 nm below the indenter– deforms as soon as the tip moves down (though the deformation initially preserves interatomic bond distances, more below). After load, the membranes are equilibrated at 1 Kelvin for 500,000 fs, with the tip remaining at a vertical distance .
iv.2 Membrane mechanics beyond first-order continuum elasticity: The isometric and anharmonic load regimes
The dimensionless quantity –with nm the closest distance from the geometrical center to the edge [Fig. 3(b)]– has been used as a measure of strain Xu et al. (2012); Klimov et al. (2012). It proves inaccurate for freestanding (rippled) membranes as the initial deformation is isometric (i.e., bond changes are initially unnoticeable). To see this, we show in Fig. 3(c) the height profiles versus , and in Fig. 3(d) the corresponding increase of the bond lengths. Even though the height plots show some amount of curvature, no significant bond length increase can be seen on the first two plots in Fig. 3(d). Indeed, when is 7% already, the largest bond increase, right below the tip, is equal to 1.2% (so that for this amount of load, the bond increase is not equal to , but rather to ). Thus, neglect of rippling Meyer et al. (2007); Fasolino et al. (2007) on thin-plate-based strain engineering (i.e., setting the initial configuration to be a plate) may lead to overestimating , an observation relevant to experimentalists generating strain on freestanding graphene with local probes. The largest bond length increase approaches for higher load, as the distortion below the tip becomes highly nonlinear (more below). Figure 3(d) also indicates bond length increases with radial symmetry near the geometrical center, determining the spatial profile of the pseudo-magnetic field that is generated by a spherical tip.
We perform an analysis of the elastic energy as a function of (Fig. 4(a)). We observe three distinct regimes. In the first regime, the elastic energy does not increase beyond fluctuations signified by error bars: This is the isometric regime, in which the initially rippled membrane follows the probe without necessarily increasing its elastic energy, nor producing significant mechanical strain. This regime holds for values of up to a few percent. The second regime in Fig. 4(a) is harmonic, as indicated by a quadratic dependence of elastic energy on . The harmonic regime holds for in a narrow range between 4 and 10%. For 10%, the system enters the anharmonic regime. We note that in the context of thin plates, only the harmonic and anharmonic regimes have been discussed in the past Duan and Huang (2009). The theory based on first-order continuum elasticity may not hold in the anharmonic regime.
The results shown in Fig. 4(a) can be understood by an analysis of the total and constituent elastic energies. The decomposition of the total elastic energy into torsional, stretching, and bending components is shown in Fig. 4(b) for triangular membrane subject to the load nm. The shaded area indicates the load time ; atomic relaxation follows in the remaining time. We observe that the total energy does not increase until nm, however the energy components provide a very interesting insight: The two leading energy contributions are the stretching and torsion of bonds. While the stretching contribution decreases –perhaps due to the fact that the tip suppresses some fluctuations in bond distances when pushing the membrane, we find that the torsion energy goes up by an almost equal amount. The remaining bending contribution to the elastic energy is an order of magnitude smaller. Thus, the total elastic energy remains practically constant for loads up to nm.
V Applying the formalism to graphene membranes under central load
v.1 Evaluation of sublattice symmetry
We plot the measures given by Equations (4-6) in Figure 5, in order to demonstrate their actual value. unit cells for which sublattice symmetry hold to numerical precision are told by the white color. As expected, deviations become larger in the close proximity of the mechanical extruder (located at the membrane’s geometrical center), and for increasing values of .
v.2 Evaluation of the flattening procedure
The vertical displacement in Fig. 2(d) is exaggerated. This displacement is less than 0.3% at distances 1 nm away from the extruder, as seen in Fig. 6. This value is one thousand times smaller than % employed to generate the atomic configuration, and represents the order of magnitude of the error introduced by the collapsing of the central atom into the x-y plane.
v.3 Gauge fields
v.4 Local density of states
The tight-binding Hamiltonian and the LDOS are meaningful regardless of the scale of the mechanical deformation. Here we display the LDOS with a 5 meV energy resolution on membranes with rims containing three million atoms. To avoid rescaling Guinea et al. (2010) we employed the Lanczos tight-binding method Wang et al. (2012).
We plot in Fig. 7(c) the LDOS with turned off at ten radial positions (see inset in Fig. 7(c)). For each radial position there are three curves, related by a 120 rotation. The curves are vertically offset for clarity. The gray v-shaped trendlines represent the DOS of unstrained graphene. We highlight a number of features: (i) A sharp zero LL, absent at some locations (points 2 and 8). (ii) Broad features in the LDOS, symmetric with respect to the zero level Levy et al. (2010); it is not clear if those correspond to a single LL or contain at least two broad LLs. From the energy locations for LLs and was estimated (assuming it uniform) and shown in some curves. Some curves are not symmetric under rotation (points 4, 6, 7, 9, and 10; not all three curves overlap). At locations 2 and 8 only a change in slope Choi et al. (2010); de Juan et al. (2012a) is seen. Only when the pseudo-magnetic field is uniform should one expect the LLs to be sharp and position-independent.
Relevance of the deformation potential in computing LDOS curves
In Fig. 7(d) we gradually turn on (Fig. 7(a)) at points 1 and 5. Importantly, equally-spaced peaks appear already at a screened , much like the equally-spaced peaks seen in Ref. Klimov et al. (2012) at zero magnetic field. is not negligible in our system as it creates a confining well.
To complete the study we plot in Fig. 7(e) the LDOS at points 1 and 5, now setting (by using a membrane with no strain) and using from Fig. 7(a): The plots in Fig. 7(d) can not be understood as a simple superposition of plots in Figs. 7(c) and those in Fig. 7(d): Both and are needed in computing the correct LDOS. Given the existence of dI/dV data obtained with local probes Klimov et al. (2012), it is important for theoretical works to report LDOS data, complementing their reported .
We have provided a theory for strain engineering valid beyond continuum elasticity, and strictly applicable for negligible curvature. We provide a measure to determine the extent to which mechanical distortions sublattice symmetry, in terms of changes in angles and lengths . For this we re-express the theory beyond continuum elasticity and explicitly on the atomic lattice. Using this formalism, we studied triangular rippled graphene membranes under mechanical load by a sharp tip. Gauge fields were computed from atomic displacements alone. We have found that rippled membranes will initially accommodate the extruder without increasing bond distances (graphene is a shell); neglecting this fact results in overestimated gauge fields. We also demonstrated in a simple way why no point dependent fields exist to first order. We studied the LDOS at many spatial locations. The scalar deformation potential gives rise to a number of equally-spaced peaks on the LDOS, even when partially screened.
Acknowledgements.We acknowledge computer support from HPC at Arkansas (RazorII), and XSEDE (TG-PHY090002, Blacklight, and Stampede), and exchanges with B. Uchoa, M. Vanević, L. Bellaiche, and M. A. Kuroda.
- J. C. Meyer, et al., Nature 446, 60 (2007).
- J. S. Bunch, et al., Science 315, 490 (2007).
- A. L. Kitt, et al., Phys. Rev. B 85, 115432 (2012).
- F. Guinea, et al., Nature Physics 6, 30 (2010).
- V. M. Pereira and A. H. Castro-Neto, Phys. Rev. Lett. 103, 046801 (2009).
- H. Suzuura and T. Ando, Phys. Rev. B 65, 235412 (2002).
- M. A. H. Vozmediano, et al., Phys. Rep. 496, 109 (2010).
- F. de Juan, A. Cortijo, M. A. H. Vozmediano, and A. Cano, Nature Phys. 7, 811 (2011).
- F. de Juan, M. Sturla, and M. A. H. Vozmediano, Phys. Rev. Lett. 108, 227205 (2012a).
- K. S. Kim, et al., Nature 457, 706 (2009).
- N. Levy, et al., Science 329, 544 (2010).
- F. de Juan, et al., Phys. Rev. B 76, 165409 (2007).
- S.-M. Choi, et al., Phys. Rev. B 81, 081407 (2010).
- A. H. Castro-Neto, et al., Rev. Mod. Phys. 81, 109 (2009).
- J. L. Ericksen, Math. Mech. Solids. 13, 199 (2008).
- N. Abedpour, et al., Phys. Rev. B 84, 115437 (2011).
- C. Pryor, et al., J. Appl. Phys. 83, 2548 (1998).
- A. Fasolino, et al., Nature Materials 6, 858 (2007).
- N. N. Klimov, et al., Science 336, 1557 (2012).
- P. Xu, et al., Phys. Rev. B 85, 121406(R) (2012).
- R. Zan, et al., Nanoscale 4, 3065 (2012).
- S. Plimpton, J. Comp. Phys. 117, 1 (1995), http://lammps.sandia.gov.
- W. H. Duan and C. M. Huang, Nanotechnology 20, 075702 (2009).
- L. E. Malbern, Introduction to Mechanics of a Continuum Medium. Prentice-Hall. (1977).
- Z. F. Wang, et al., Nano Lett. 12, 3833 (2012).
- C.-H. Park, et al., Nature Phys. 4, 213 (2008).