Pseudomagnetic fields in graphene nanobubbles of constrained geometry: A molecular dynamics study
Zenan Qi, Alexander L. Kitt, Harold S. Park, Vitor M. Pereira, David K. Campbell, A. H. Castro Neto
PPseudomagnetic fields in graphene nanobubbles of constrained geometry: A moleculardynamics study
Zenan Qi ∗ , Alexander L. Kitt † , Harold S. Park ‡ , § Vitor M. Pereira ¶ , ∗∗ David K. Campbell †† , and A. H. Castro Neto ‡‡
2, 3, 4 Department of Mechanical Engineering, Boston University, Boston, MA 02215 Department of Physics, Boston University, 590 Commonwealth Ave, Boston, Massachusetts 02215, USA Graphene Research Centre & Department of Physics,National University of Singapore, 2 Science Drive 3, Singapore 117542 Department of Electrical and Computer Engineering,National University of Singapore, 4 Engineering Drive 3, Singapore 117583 (Dated: August 13, 2018)Analysis of the strain-induced pseudomagnetic fields generated in graphene nanobulges underthree different substrate scenarios shows that, in addition to the shape, the graphene-substrate in-teraction can crucially determine the overall distribution and magnitude of strain and those fields,in and outside the bulge region. We utilize a combination of classical molecular dynamics, contin-uum mechanics, and tight-binding electronic structure calculations as an unbiased means of study-ing pressure-induced deformations and the resulting pseudomagnetic field distribution in graphenenanobubbles of various geometries. The geometry is defined by inflating graphene against a rigidaperture of a specified shape in the substrate. The interplay among substrate aperture geometry,lattice orientation, internal gas pressure, and substrate type is analyzed in view of the prospect ofusing strain-engineered graphene nanostructures capable of confining and/or guiding electrons atlow energies. Except in highly anisotropic geometries, the magnitude of the pseudomagnetic field isgenerally significant only near the boundaries of the aperture and rapidly decays towards the centerof the bubble because under gas pressure at the scales considered here there is considerable bend-ing at the edges and the central region of the nanobubble displays nearly isotropic strain. Whenthe deflection conditions lead to sharp bends at the edges of the bubble, curvature and the tiltingof the p z orbitals cannot be ignored and contributes substantially to the total field. The strongand localized nature of the pseudomagnetic field at the boundaries and its polarity-changing pro-file can be exploited as a means of trapping electrons inside the bubble region or of guiding themin channel-like geometries defined by nano-blister edges. However, we establish that slippage ofgraphene against the substrate is an important factor in determining the degree of concentrationof pseudomagnetic fields in or around the bulge since it can lead to considerable softening of thestrain gradients there. The nature of the substrate emerges thus as a decisive factor determiningthe effectiveness of nanoscale pseudomagnetic field tailoring in graphene. PACS numbers: 81.05.ue, 73.22.Pr, 71.15.Pd, 61.48.Gh
I. INTRODUCTION
Since the discovery of a facile method for its isolation,graphene , the simplest two-dimensional crystal, has at-tracted intense attention not only for its unusual physicalproperties , but also for its potential as the basic build-ing block for a wealth of device applications. There existkey limitations that appear to restrict the applicationof graphene for all-carbon electronic circuits: one suchlimitation is that graphene, in its pristine form, is wellknown to be a semi-metal with no band gap . A highlyactive field of study has recently emerged based on theidea of applying mechanical strain to modify the intrinsic ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] ¶ Electronic address: [email protected] †† Electronic address: [email protected] ‡‡ Electronic address: [email protected] response of electrons to external fields in graphene .This includes the strain-induced generation of spectral(band) gaps and transport gaps, which suppress conduc-tion at small densities. In this context, several groups have employed continuum mechanics coupled witheffective models of the electronic dynamics to study thegeneration of pseudomagnetic fields (PMFs) in differentgraphene geometries and subject to different deforma-tions. The potential impact of strain engineering beyondthe generation of bandgaps has also attracted tremen-dous interest .Pereira et al. showed that a band gap will not emergeunder simple uniaxial strain unless the strain is largerthan roughly 20 %. This theoretical prediction, based onan effective tight-binding model for the electronic struc-ture, has been subsequently confirmed by various moreelaborate ab-initio calculations . The robustness ofthe gapless state arises because simple deformations ofthe lattice lead only to local changes of the position ofthe Dirac point with respect to the undeformed latticeconfiguration and to anisotropies in the Fermi sur- a r X i v : . [ c ond - m a t . m e s - h a ll ] S e p face and Fermi velocity . The shift in the position ofthe Dirac point is captured, in the low-energy, two-valley,Dirac approximation, by a so-called pseudomagnetic vec-tor potential and resulting pseudomagnetic field (PMF)that arises from the strain-induced perturbation of thetight-binding hoppings . As a result, electrons reactto mechanical deformations in a way that is analogousto their behavior under a real external magnetic field,except that overall time-reversal symmetry is preserved,since the PMF has opposite signs in the two time-reversalrelated valleys .Guinea et al. found that nearly homogeneous PMFscould be generated in graphene through triaxial stretch-ing, but the resulting fields were found to be moder-ate, unless relatively large ( i.e. , >
10 %) tensile strainscould be applied. Unfortunately, such large planar tensilestrains have not been experimentally realized in grapheneto date. This is arguably attributed to the record-hightensile modulus of graphene and the unavoidable diffi-culty in effectively transferring the required stresses fromsubstrates to this monolayer crystal .It is thus remarkable that recent experiments re-port the detection of non-uniform strain distributions inbubble-like corrugations that generate PMFs locally ho-mogeneous enough to allow the observation of Landauquantization by local tunneling spectroscopy. The mag-nitude of the PMFs reported from the measured Lan-dau level spectrum reaches hundreds (300 to 600) of Tes-las , providing a striking glimpse of the impact thatlocal strain can potentially have on the electronic prop-erties. A difficulty with these experiments is that, upto now, such structures have been seen and/or gener-ated only in contact with the metallic substrates thatare used in the synthesis of the sample. This is an obsta-cle, for example, to transport measurements, since thiswould require the transfer of the graphene sheet to an-other substrate, thereby destroying the favorable localstrain distribution. In addition, a systematic study of dif-ferent graphene bubble geometries and substrate types,which could reveal the subtleties that different geome-tries bring to the related strain-induced PMFs has notbeen reported. Furthermore, most previous studies ofthe interplay between strain and electronic structure ingraphene have addressed the deformation problem froman analytic continuum mechanics point of view, with theexception of a few recent computational studies .It is in this context that we report here resultsfrom classical molecular dynamics (MD) simulations ofstrained graphene nanobubbles induced by gas pres-sure. The MD simulations are used to complement andcompare continuum mechanics approaches to calculatingstrain, in order to examine the pressure-induced PMFsin ultra-small graphene nanobubbles of diameters on theorder of 5 nm. Controlled synthesis of such small strainednanobubbles has gained impetus following the recent ex-periments by Lu et al. . Our aim is to use an unbi-ased calculation for the mechanical response of grapheneat the atomistic level, on the basis of which we can (i) FIG. 1: (Color online) Illustration of the strategy employed inour studies to generate nanobubbles by pressurizing graphenethrough a predefined substrate aperture. The picture showsone of the actual simulation cells used in our MD compu-tations. In gold, gray and red colors are represented, respec-tively, the Au substrate, the graphene sheet and the Ar atoms.A hole is carved in the Au substrate (perimeter outlined), andits perimeter geometry determines the shape of the resultinggraphene bubble. Visualization is performed using VMD . extract the relaxed lattice configurations without any as-sumptions; (ii) calculate the PMF distribution associatedwith different nanobubble geometries; (iii) discuss the in-fluence of substrate and aperture shape on PMF distri-bution; (iv) identify conditions under which explicit con-sideration of the curvature is needed for a proper accountof the PMFs.We first describe the simulation methodology that wasemployed to determine the atomic displacements fromwhich the strain tensor, modified electronic hopping am-plitudes, and PMFs can be obtained. This is followed bynumerical results of the strain-induced PMFs for differ-ent graphene nanobubble geometries in a simply clampedscenario. We next discuss the considerable importance ofthe substrate interaction and, finally, analyze the relativecontributions of orbital bending and bond stretching tothe total PMF. II. SIMULATION METHODOLOGY
Recent experiments have shown that graphenenanobubbles smaller than 10 nm can be prepared onmetallic substrates, and that large PMFs in the hundredsof Tesla result from the locally induced non-homogeneousstrain . Because such small nanobubbles can be di-rectly studied using classical MD simulations, we employMD to obtain the deformed graphene bubble configura-tions due to an externally applied pressure.The atomisticpotentials that describe the carbon-carbon interactionshave been extensively investigated and, hence, graphene’snano-mechanics can be simulated without any particularbias, and to a large accuracy within MD. Once the de-formation field is known from the simulations, we obtainthe strain distribution in the inflated nanobubble, finallyfollowed by a continuum gauge field approach to extractthe resulting PMF distribution . A. Details of the MD simulations
Our MD simulations were done with the Sandia-developed open source code LAMMPS . Thegraphene nanobubble system consisted of three parts, asillustrated in Fig. 1: a graphene monolayer, the hexag-onal (111) surface of an FCC gold substrate, and ar-gon gas which was used to inflate the graphene bub-ble. We used the AIREBO potential to describe theC-C interactions, as this potential has been shown to de-scribe accurately the various carbon interactions includ-ing bond breaking and reforming . The substrate-graphene and gas-graphene interactions were modeled bya standard 12-6 Lennard-Jones potential: V ( r ij ) = 4 (cid:15) ij (cid:34)(cid:18) σ ij r ij (cid:19) − (cid:18) σ ij r ij (cid:19) (cid:35) , (1)where r ij represents the distance between the i -th carbonand the j -th gold atom.The dimension of the simulation box was20 × × , and the substrate was comprised ofAu atoms with a thickness of 2 nm, or about 2.5 timesthe cutoff distance of the interatomic potential . Aper-tures of different shapes (viz. triangle, rectangle, square,pentagon, hexagon, and circle) were “etched” in thecenter of the substrate to allow the graphene membraneto bulge inwards due to the pressure exerted by the Argas. The whole system was first relaxed for 50 ps, atwhich time the Ar gas was pushed downward (as in apiston) to exert pressure on the graphene monolayer,causing it to bulge inward in the shape cut-out from thegold substrate. The system is then allowed to equilibrateagain under the increased gas pressure. All simulationswere carried out at room temperature (300 K) usingthe Nose-Hoover thermostat . The choice of Ar in ourcalculations is not mandatory. Substitution with othermolecular species should pose no difficulty, the samebeing true regarding the substrate, as shown previouslyin references 41 and 34.To elucidate the effect of different substrates on thePMF distributions in the nanobubbles, we perform MDsimulations with two different substrates, in addition toperforming the simulations with fixed edges and no sub-strate. Specifically, we used both Au and Cu (111) sub-strates, where the detailed parameters and descriptionswill be discussed in later sections.After obtaining the graphene bubble, we held thepressure constant for 10 ps to achieve thermal equilib-rium. We note that during the entire simulation nogas molecules leaked away from the system, which againdemonstrates the experimentally observed atomic imper-meability of monolayer graphene . Our simulations are close in spirit to the experimentsreported in reference 45, but targeting smaller hole aper-tures due to computation limitation. We note that thismethod of using gas-pressure to generate the graphenenanobubbles is different from the situations explored inthe recent experiments that focus on the PMF distri-bution . However, it is in some ways more control-lable due to the utilization of a substrate with a distinctpattern coupled with externally applied pressure to forcegraphene through the patterned substrate to form a bub-ble with controllable shape and height.The final (inflated bubble) configuration gives us thebasic ingredients needed to extract the strain distributionin the system, as well as the perturbed electronic hop-ping amplitudes. To calculate the strain directly fromthe displaced atomic positions we employ what we shalldesignate as the displacement approach . We note thata previous study used a stress approach for a similarcalculation. However, the stress approach fails to predictreasonable results in our case, which we attribute to theinability of the virial stresses to properly convey the to-tal stress at each atom of the graphene sheet when theload results from interaction with gas molecules. Further-more, in the stress approach one assumes a planar (and,in addition, usually linear) stress-strain constitutive rela-tion which leads to errors when large out-of-plane defor-mations arise, as in the case of the nanobubbles. Furtherdetails on the strain calculation are given in appendix D. B. Displacement approach to calculate strain
In continuum mechanics the infinitesimal strain tensoris written in Cartesian material coordinates ( X i ) as (cid:15) ij = 12 (cid:18) ∂u i ∂X j + ∂u j ∂X i (cid:19) + 12 (cid:18) ∂u k ∂X i ∂u k ∂X j (cid:19) . (2)To utilize Eq. 2, it is clear that the displacement fieldmust be obtained such that its derivative can be eval-uated to obtain the strain. In order to form a linearinterpolation scheme using finite elements , we exploitthe geometry of the lattice and mesh the results of ourMD simulation of the deformed graphene bubble usingtetrahedral finite elements defined by the positions offour atoms: the atom of interest (with undeformed co-ordinates R ), and its three neighbors (with undeformedcoordinates R , R , R ). After deformation, the newpositions of the atoms are r , r , r and r , respectively.To remove spurious rigid body translation and rotationmodes, we took the atom of interest ( R ) as the refer-ence position, i.e. , r = R . The displacement of its threeneighbors could then be calculated, and subsequently thecomponents of the strain tensor (cid:15) ij were obtained by nu-merically evaluating the derivative of the displacementinside the element. C. Pseudomagnetic fields
Non-zero PMFs arise from the non-uniform displace-ment in the inflated state. These PMF reflect the phys-ical perturbation that the electrons near the Fermi en-ergy in graphene feel as a result of the local changes inbond length. It emerges straightforwardly in the follow-ing manner. Nearly all low-energy electronic propertiesand phenomenology of graphene are captured by a sim-ple single orbital nearest-neighbor tight-binding (TB) de-scription of the π bands in graphene . In second quan-tized form this tight-binding Hamiltonian reads H = − (cid:88) i, n t (cid:0) r i , r i + n (cid:1) a † r i b r i + n + H. c. , (3)where t (cid:0) r i , r i + n (cid:1) represents the hopping integral be-tween two neighboring π orbitals, n runs over the threenearest unit cells, and a r i ( b r i ) are the destruction op-erators at the unit cell r i and sublattice A ( B ). In theundeformed lattice the hopping integral is a constant: t (cid:0) r i , r i + n (cid:1) = t (cid:0) R i , R i + n (cid:1) = t = 2 . t ij between neighboring atoms i and j , which is exponen-tially sensitive to the interatomic distance. The othereffect is caused by the curvature induced by the out-of-plane deflection, which means that the hopping ampli-tude is no longer a purely V ppπ overlap (in Slater-Kosternotation) but a mixture of V ppπ and V ppσ . More pre-cisely, one can straightforwardly show that the hoppingbetween two p z orbitals oriented along the unit vectors n i and n j and a distance d apart is given by − t ij = V ppπ ( d ) n i · n j + V ppσ ( d ) − V ppπ ( d ) d ( n i · d )( n j · d ) . (4)To capture the exponential sensitivity of the overlap in-tegrals to the interatomic distance d we model them by V ppπ ( d ) = − t e − β ( d/a − , (5a) V ppσ ( d ) = +1 . t e − β ( d/a − , (5b)with a (cid:39) .
42 ˚A the equilibrium bond length in graphene.For static deformations a value β ≈ V ppπ ( d ) in agreement withfirst-principles calculations ; we use the same decayconstant β for both overlaps, which is justified from aM¨ulliken perspective since the principal quantum num-bers of the orbitals involved is the same .In the undeformed state Eq. (4) reduces to − t ij = V ppπ ( a ) ≡ − t and is, of course, constant in the entiresystem. But local lattice deformations cause t (cid:0) r i , r i + n (cid:1) to fluctuate, which we can describe by suggestivelywriting t (cid:0) r i , r i + n (cid:1) = t + δt (cid:0) r i , r i + n (cid:1) . In the low energy(Dirac) approximation, the effective Hamiltonian aroundthe point ± K in the Brillouin zone can then be writtenas H ± K eff = v F σ · (cid:0) p ∓ q A (cid:1) , (6)where (cid:126) v F = 3 ta/ q represents the charge of the currentcarriers ( q > q < A = A x e x + A y e y are given explicitly in termsof the hopping perturbation by A x ( R ) − iA y ( R ) = 1 qv F (cid:88) n δt (cid:0) r , r + n (cid:1) e i K · n . (7)For nearly planar deformations (small out-of-plane vs in-plane displacement ratios and thus neglecting bendingeffects) δt can be expanded in terms of the local dis-placement field and, consequently, can be cast in termsof the strain components. Orienting the lattice so thatthe zig-zag direction is parallel to e x leads to A x ( R ) − iA y ( R ) (cid:39) (cid:126) β qa (cid:0) (cid:15) xx − (cid:15) yy + 2 i (cid:15) xy (cid:1) , (8)Since we are ultimately interested in the PMF, only thecontributions to A ( R ) arising from the hopping modifi-cation are considered here, as they are the ones that sur-vive after the curl operation ; we also don’t con-sider contributions beyond second order smallness ( ∼ k (cid:15) , ∼ k , etc.). In the planar strain situation the whole in-formation about the electronic structure is reduced to theparameter β = − ∂ log t ( r ) /∂ log r (cid:12)(cid:12) r = a .From the coupling in Eq. (6) where the effects of strainare captured by replacing p → p − q A it is clear that thelocal strain is felt by the electrons in the K valley inthe same way as an external magnetic field would be.In particular, we can quantify this effect in terms of thePMF, which is defined as B = ∂ x A y ( R ) − ∂ y A x ( R ) . (9)This is the central quantity of interest in this work; inthe subsequent sections the combined effects of gas pres-sure, hole geometry, and substrate interaction will be an-alyzed from the point of view of the resulting magnitudeand space distribution of the PMF, B , obtained in thisway. For definiteness we set q = e , e being the elemen-tary charge, which means that we shall be analyzing thePMF from the perspective of holes ( q > B can be calculated directly fromEq. (7) by computing the hopping between all pairs ofneighboring atoms in the deformed state, or from Eq. (8)by calculating the strain components throughout the en-tire system as described in the previous section. Theformer strategy is here referred to as the TB approach ,and the latter as the displacement approach , as per thedefinitions in section II B. Our PMF calculations in thefollowing sections are done by following the TB approach,except when we want to explicitly compare the results ob-tained with the two approaches. In those cases, such as inthe next section or in Appendix C, that will be explicitlystated.
III. CLAMPED GRAPHENE NANOBUBBLES
We first simulated an idealized system consisting onlyof Ar gas molecules and graphene, neglecting the interac-tion with the underlying substrate, and where we strictlyfixed all carbon atoms outside the aperture region duringsimulation. This provides a good starting point to un-derstand how the shape of the substrate aperture affectsthe PMF distribution. A similar system has been usedin previous work , as this corresponds to a continuummodel with clamped edges .We start with the most symmetric geometry, a circulargraphene bubble, and compare the atomistic result withthe continuum Hencky solution . In contrast to smalldeformation continuum models , the Hencky model isvalid for large in-plane (stretching) deformations, whichlead to a different PMF distribution. To compute thePMFs associated with this analytical solution we usedEq. (8). Figs. 2(c,d) show that the PMF distributionis dominated by very large magnitudes at the edges fol-lowed by a rapid decay towards the inside region of thenanobubble. Both the MD and Hencky results show thesix-fold symmetry expected for a cylindrically symmet-ric strain distribution; this agreement demonstrates theMD simulation successfully captures the strain distribu-tion underlining the computed PMF. There are, however,two quite clear discrepancies between the PMF in thesetwo figures: (i) Hencky’s solution (panel d) yields valuesconsiderably smaller in magnitude than the calculationbased on the MD deformations combined with Eqs. (4)and (7) (panel c); (ii) the sign of the PMF in panel (d) is apparently reversed with respect to the sign of panel (c).These discrepancies stem from the substantial bendingpresent in graphene near the hole perimeter, and deservea more detailed inspection in terms of the relative mag-nitude of the two contributions to the hopping variation:bond stretching and bond bending.Since Hencky’s result of Fig. 2(d) hinges on Eq. (8)that expresses the vector potential directly in terms ofthe strain tensor components, let us start by analyzingthe predictions obtained by applying it to the atomisticcase as well; to do that one computes the strain fromthe MD simulations using the displacement approach dis-cussed earlier. The result of that is shown in Fig. 2(f),where the most important difference in comparison withFig. 2(c) is the significant reduction of the maximal fieldsobtained near and at the edges; this reflects the error in-curred in the quantitative estimate of B when the effect ofbending is neglected. Note that, by construction, Eq. (8)accounts only for the bond-stretching, and is accurateonly to linear order in strain because it is based on a lin- ear expansion of the hopping in the interatomic distances.Hence, in order to correctly extract from the atomisticsimulations the total stretching contribution beyond lin-ear order while still ignoring bending effects, we shouldcalculate the PMF with the hopping as defined in Eq. (4)(TB approach), but explicitly setting n i · n j = 1 and n i · d = 0 ( i.e. assuming local flatness). The outcome ofthis calculation is shown in Fig. 2(e) which, in practi-cal terms, is the counterpart of Fig. 2(c) with bendingeffects artificially suppressed. In comparison with panel(f), it leads to slightly smaller PMF magnitudes. Thelinear expansion in strain of Eq. (7) thus slightly overesti-mates the field magnitudes, something expected becausethe hopping is exponentially sensitive to the interatomicdistance and, by expanding linearly, one overestimatesits rate of change with distance, overestimating the fieldmagnitude as a result. One key message from Fig. 2 andthe comparison between panel (c) and any of the subse-quent ones is that the effects of curvature are significantat these scales of deflection and bubble size, particularlyat the edge, where they clearly overwhelm the “in-plane”stretching contribution. We will revisit this in more de-tail in section VI.The second key message gleaned from Fig. 2 pertains tothe importance of properly considering the boundary andloading conditions when analytically modeling the strainand deflection of graphene. This is related to the appar-ent opposite sign in the PMF at the edge obtained fromHencky’s solution in panel (d) when compared with allthe other panels (containing the MD-derived results). Toelucidate the origin of the difference we show in Fig. 3(a)the PMF divided by the angular factor sin(3 θ ), and aver-aged over all the angles (details discussed in appendix A).This plot provides a summary of the data in Figs. 2(d,e,f)and allows a cross-sectional view of the variation of thefield magnitude with distance from the center of thenanobubble. Direct inspection shows that the averagedMD data follows Hencky’s prediction inside the bubblenearly all the way to the edge, at which point the PMFderived from the atomistic simulations swerves sharplyupwards, changes sign, and returns rapidly to zero withinone lattice spacing beyond the bubble edge (the curvederived from Hencky’s model terminates at the edge, byconstruction). This effective sectional view explains whythe density plots in Figs. 2(c,d) seem to have an overallsign mismatch: in the MD-derived data, the plots of thePMF distribution are dominated by the large values atthe edge which have an opposite sign to the field in theinner region. Fig. 3(a) shows that, rather than a discrep-ancy, there is a very good agreement between the strainfield predicted by Hencky’s solution and a fully atom-istic simulation throughout most of the inner region ofthe nanobubble. However, since Hencky’s solution as-sumes fixed boundary conditions at the edge (zero de-flection, zero bending moment) , it cannot capture thesharp bends expected at the atomic scale generated bythe clamping imposed in these particular MD simulations(in effect, corresponding to zero deflection and its deriva- FIG. 2: (Color online) Results for a circular graphene bubble with 4 nm radius and pressurized up to ∼ with the axes scaledin units of the circle radius, (e) PMF by TB method with in-plane component only, (f) PMF by displacement method. Notethat, except for (d), all the panels refer to the same atomistic configuration. PMF in shown in units of Tesla. The edge of thesubstrate aperture used in the MD simulation is outlined (gray line) for reference. tive). The finite bending stiffness of graphene comesinto play in that region, generating additional strain gra-dients which explain the profile and large magnitude ofthe PMF seen in the atomistic simulations.In Fig. 3(b) we plot the evolution of the deflection andmaximum PMF with increasing gas pressure. The max-imum PMF is obtained around the edge of the aperture,and the values shown in the figure correspond to an angu-lar average of the PMF amplitude there (see appendix Afor details). The MD and analytical (Hencky’s) solutionsgive comparable results for the deflection in the pressurerange below < × bar (Fig. 3(b), right vertical scale).At higher pressures, Figs 3(b) and 3(c) show that the an-alytical solution yields a slightly smaller deflection, as theunderlying model does not capture the nonlinear elasticsoftening that has been observed in graphene in both ex-periments and previous MD simulations . Fig. 3(b)includes also the maximum PMFs occurring at the bub-ble edge, when computed with the different approachesdiscussed above in connection with Figs 2(c-f). We high-light that Hencky’s solution cannot generate significantPMFs even at the largest deflections, whereas experi-ments in similarly sized and deflected nanobubbles easilyreveal PMFs in the hundreds of Teslas . This raises questions about the applicability of the Hencky solutionat these small scales and large deflections.The pressure required to rupture this graphene bub-ble was determined to be around 1 . × bar from ourMD simulations. Such a large value is required becauseof the small dimensions of the bubble. We can calculatethe fracture stress by adopting a simple model for a cir-cular bulge test, i.e. , σ ∼ RδP w , where σ , δP , R , and w are the stress, pressure difference, radius, and thicknessof the membrane, respectively. Assuming w to be 3.42 ˚A,we obtain a fracture strength of about 80 GPa, which isin agreement with previous theoretical and experimen-tal results. Note that the plot in Fig. 3 shows verylarge pressures (up to near the rupture limit of the bub-ble) and correspondingly large deflections since we wishto highlight the points of departure between the elasticmodel and the simulation results. Pressures and deflec-tions considered in the specific cases discussed below areconsiderably smaller.With the good performance of the atomistic model onthe circular graphene bubble established, we next extendthe analysis to nanobubbles with different shapes. Thebubbles are similarly obtained by inflation of grapheneunder gas pressure against a target hole in the substrate (cid:72) Å (cid:76) P M F (cid:144) S i n (cid:72) Θ (cid:76)(cid:72) T (cid:76) DispTB, in (cid:45) planeHencky (a)
PMF/Sin(3 q ) (T) P r e s s u r e ( k B a r ) P M F ( D i s p ) P M F ( T B , i n - p l a n e ) P M F ( T B , c o m p l e t e , x 0 . 1 ) P M F ( H e n c k y )
D e f l e c t i o n ( M D ) D e f l e c t i o n ( H e n c k y )
Deflection ( (cid:1) ) (b) −40 −20 0 20 400246810 Radius (Å) D e f l e c t i on ( Å ) Diametral Height
MDHencky (c)
FIG. 3: (Color online) (a) Angular-averaged values of B/ sin(3 θ ) for the circular nanobubble with R = 4 nm con-sidered in Fig. 2. The different datasets correspond to differ-ent strategies discussed in the main text to obtain the PMF.The vertical line at r = R ≈
40 ˚A marks the radius of thecircular aperture in the substrate. For r < R the results ex-tracted from MD closely follow the analytical curve, but thereis a sharp sign change and increase at r ≈ R (see appendixA for details of the averaging procedure, as well as for theTB data including the full hopping perturbation). (b) Com-parison between the pressure-induced deflection and maxi-mum PMF magnitude at the edge, | B ( R ≈
40 ˚A) | , obtainedwith the different approximations discussed in the text. Thepoints corresponding to the complete TB hopping are scaledby 0.1 for better visualization. (c) A section of the simu-lated nanobubble (MD) at ∼
19 kBar and the correspondingHencky’s solution (the inset shows a 3D perspective of the for-mer with the color scale reflecting the vertical displacement). with the desired shape. Fig. 4 shows results of a studyof different shapes to which the displacement interpola-tion approach was applied to obtain the strain field and,thus, the PMFs. The shapes are a square, a rectangle(aspect ratio of 1:2), a pentagon, a hexagon, and a circle,and are presented in order of approximately decreasingsymmetry. Those geometries are chosen because they aresufficiently simple that they can be readily fabricated ex-perimentally with conventional etching techniques. Thedimensions of the different bubbles were chosen such thattheir areas were approximately ∼
50 nm . The pressurewas 19000 bar and side lengths for the bubble geometriesshown in Figs. 4(a)-4(f) were, respectively, 4 nm (circle),4.4 nm (hexagon), 5.7 nm (pentagon), 5 nm (rectangle,short edge), 7.1 nm (square), 10.6 nm (triangle).It is worth emphasizing that these features depend onthe orientation of the graphene lattice with respect to thesubstrate aperture, as we would expect. This is clearlyvisible in the case of the square bubble in Fig. 4(e), forwhich the sharp magnetic field along the boundary ispresent along the horizontal (zig-zag) edges of the bub-ble but not along the vertical ones (armchair). This isalso the reason why only the triangular aperture shown inFig. 4(f) leads to a strong PMF that is nearly uniform asone goes around the boundary of the nanobubble. This isan important consideration for the prospect of engineer-ing strained graphene nanostructures capable of guidingor confining electrons within, much like a quantum dot .The sharp PMF at the boundary acts effectively as astrong magnetic barrier, which might be tailored to con-fine some of the low energy electronic states .The resulting PMF patterns in Fig. 4 show that thehighest values are found at the corners and edges of thedifferent bubble shapes. To illustrate more clearly thePMF patterns, we inflated the bubbles to large deflec-tions ( ∼ × bar. These largedeflections explain why the PMF magnitudes in Fig. 4may reach over 500 T. Given that the gas pressures usedto achieve the results shown in this figure are rather high,some comments are in order.First, we emphasize that the relevant parameter is thedeflection, rather than the pressure itself. In other words,gas pressure was employed here as one way of generat-ing graphene nanobubbles with predefined boundary ge-ometries and target deflections, but other loading condi-tions might be used to achieve the same parameters. Ourchoice is motivated by the desire to constrain grapheneand its interaction with the substrate as little as possi-ble. Since we intend to reproduce bubbles with lateralsize and deflections matching the magnitude of the valuesobserved experimentally this requires large pressures(for a given target deflection P is naturally smaller forlarger apertures). Secondly, Lu et al. reported thatexperimental bond elongations, estimated from directSTM mapping of the atomic positions and deflections,can exceed 10 % in graphene nanobubbles on Ru. Thehigh pressures considered in our MD simulations allow (a) (b) (c) (d) (e) (f) FIG. 4: (Color online) Top views of PMF patterns for graphene bubbles of different geometries without substrate. (a) circle(b) hexagon (c) pentagon (d) rectangle (aspect ratio 1:2) (e) square (f) triangle. All the bubble areas are ∼
50 nm , and sidelengths and pressures can be found in the main text. In all cases, the graphene lattice is oriented with the zig-zag directionalong the horizontal. The same color scale (in Tesla) is used in all panels. The edge of the substrate apertures used in the MDsimulations is outlined (gray line) for reference. us to reach bond elongations of this order of magnitude.Thirdly, pressures of the order of 10 kbar (1 GPa) havebeen recently estimated to occur within nanobubbles ofsimilar dimensions and deflections to the ones consid-ered here, formed upon annealing of graphene-diamondinterfaces . Thus, pressures of this magnitude are notunrealistic in the context of nanoscale graphene blisters. IV. SUBSTRATE INTERACTION: GRAPHENEON AU (111)
Having considered the ideal case of graphene with-out a substrate, we move forward to study the morerealistic case of graphene lying on an Au (111) sub-strate. The main difference is that the carbon atomsare not rigidly attached to the substrate anymore out-side the aperture, meaning that graphene can slideinto the aperture during inflation, subject to the in-teraction with the substrate. This is an importantqualitative difference, and reflects more closely theexperimental situation, as recently reported in refer-ence 61. The interatomic interactions were param-eterized with (cid:15) C − Au =0.02936 eV, σ C − Au =2.9943 ˚A ; (cid:15) C − Ar =0.0123 eV, σ C − Ar =3.573 ˚A ; (cid:15) Ar − Ar =0.0123 eV, σ Ar − Ar =3.573 ˚A ; the Ar-Au (gas-substrate) interac-tions were neglected to save computational time, and thesubstrate layer was held fixed for the entire simulationprocess. Most of the graphene layer was unconstrained,except for a 0.5 nm region around the outer edges of thesimulation box where it remained pinned. Since the in-teraction with the substrate is explicitly taken into con-sideration, this approach realistically describes the slid-ing and sticking of graphene on the substrate as the gaspressure is increased, as well as details of the interactionwith the substrate in and near the hole perimeter.We start the discussion with a direct comparison of thedeformation state of a circular bubble obtained from oursimulations with the predictions of a recently developedand experimentally verified ‘extended-Hencky’ model that accounts for the same sliding and friction effects.As can be seen in Fig. 5(a), after fitting the frictionin the continuum model to the MD simulation there isa very good agreement between the MD and extendedHencky results for the radial and tangential strains, (cid:15) rr and (cid:15) θθ , both in the inner and outer regions with respectto the substrate aperture. The same good agreementis seen in the PMF profile extracted from the MD andanalytical approaches, which is presented in Fig. 5(b).The numerical data points shown in this panel represent FIG. 5: (Color online) (a) Strain components (cid:15) rr and (cid:15) θθ of agraphene bubble pressurized to a deflection of ∼ (solid line) and from MD simulations within the TB(blue) or displacement (red) approach. Panel (d) shows theangular dependence of the PMF for selected radii. an angular average over an annulus centered at differ-ent radii. An important message from Fig. 5(b) is thatthe maximum magnitude of the PMF occurs around theedge of the aperture, but on the outside of the bubble.Whereas one expects the maximal PMFs to occur aroundthe edge where the strain gradients are larger, the factthat the magnitude is considerably higher right outsiderather than inside is not so obvious. This has importantimplications for the study of PMFs in graphene nanos-tructures but has been ignored by previous studies. Itimplies that models where only the deflection inside theaperture is considered (such as the simple Hencky model)can miss important quantitative and qualitative features.They are captured here because the friction and slidingeffects due to graphene-substrate interactions are natu-rally taken into account from the outset. One conse-quence is the “leakage” of strain outside the bubble re-gion and the concurrent emergence of PMFs outside theaperture. This should be an important consideration indesigning nanoscale graphene devices with functionalitiesthat rely on the local strain or PMF distribution.The other shapes studied on the Au (111) substrate areshown in Fig. 6. The dimensions are the same as in Fig. 4,with an applied pressure of ∼
30 kbar. In addition to theappearance of non-negligible PMF outside the apertureregion, a comparison with the data for bubbles clampedto the hole perimeter shows that now the PMF distri-bution inside is noticeably perturbed, and that the largefield magnitudes observed in Fig. 4 along the perimeterare considerably reduced and smoother.To understand the origin of this difference, let us an-alyze in detail the representative case of a triangularnanobubble, as previous experiments have shown thatsuch nanobubbles can exhibit PMFs in excess of 300 T . Using our MD-based simulation approach, we calculatethe PMFs for triangular graphene bubbles by inflatinga graphene monolayer through a triangular hole in thesubstrate. The set-up is as illustrated in Fig. 1, but withthe circular hole replaced by a triangular one. The trian-gular hole in the substrate had a side length of 10.6 nm,and the graphene sheet was inflated to a deflection of ∼ , this tri-axial symmetry is crucialfor the experimental observation of Landau levels in ref-erence 31 because it leads to a quasi-uniform PMF insidethe nanobubble. Inspection of Fig. 4(f) confirms thatthe field is indeed of significant magnitude and roughlyuniform within the bubble.When the full interaction with the substrate is includedand the graphene sheet is allowed to slip and slide to-wards the aperture under the inflation pressure, the ge-ometry is no longer as effective as before in generatinga clear triaxial symmetry: a comparison of the top andbottom rows of Fig. 7 shows that the triaxial symme-try of the strain distribution is not so sharply definedin this case. Therefore, the finite and roughly uniformPMF inside the triangular boundary that is seen clearlyin Figs. 4(f) [and also Figs. 16(f)] is largely lost here. Tounderstand the difference we start by pointing out thatthe orientation of the triangular hole with respect to thecrystallographic axes used here is already the optimumorientation in terms of PMF magnitude, with its edgesperpendicular to the (cid:104) (cid:105) directions ( i.e. , parallel to thezig-zag directions). Secondly, since the graphene sheetis allowed to slide, the strain distribution in the centralregion of the inflated bubble tends to be more isotropic,as we expect for an inflated membrane because of theout of plane displacement, and as can be seen in Fig. 7.This means that the trigonal symmetry imposed on theoverall strain distribution by the boundaries of the holeis less pronounced near the center. As a result, eventhough strain increases as one moves from the edge to-wards the center (as measured, for example, by looking atthe bond elongation directly from our MD simulations),the magnitude of the PMF decreases because the trigonalsymmetry and strain gradients become increasingly lesspronounced, and we know that the isotropic (circular)hole yields zero PMF at the apex (Fig. 2).The differences in trend and the sensitivity of the PMFdistribution to the details of the interaction with the sub-strate highlight the importance of the latter in determin-ing the final distribution and magnitude of the PMF, inaddition to the loading, hole shape, and boundary condi-tions. In order to stress this aspect, and to make the roleof the substrate interaction even more evident, we shallconsider next a different metal surface.0 (a) (b) (c) (d) (e) (f) FIG. 6: (Color online) Top views of PMF patterns for graphene bubbles of different geometries on Au (111) substrates. (a)circle (b) hexagon (c) pentagon (d) rectangle (aspect ratio 1:2) (e) square (f) triangle. All the bubble areas are ∼
50 nm , andside lengths and pressures can be found in the main text. In all cases, the graphene lattice is oriented with the zig-zag directionalong the horizontal. The same color scale (in Tesla) is used in all panels. The edge of the substrate apertures used in the MDsimulations is outlined (gray line) for reference.FIG. 7: (Color online) Spatial patterns of the strain tensorcomponents (cid:15) rr and (cid:15) θθ for a triangular bubble with a 10.6 nmside. (a) and (b) pertain to graphene on a Au (111) substratewhose PMF profile has been shown in Fig. 4(f), while (c) and(d) correspond to the graphene bubble with an artificiallyfixed boundary condition whose PMF is shown in Fig. 6(f).The edge of the substrate aperture used in the MD simulationis outlined (gray line) for reference. V. SUBSTRATE INTERACTION: GRAPHENEON CU (111)
To gain further insight into the important effects ofsubstrate interactions, we carried out simulations for aCu (111) substrate, in addition to the Au (111) case con-sidered above. This is in part motivated by a recentexperimental study showing that graphene grown bychemical vapor deposition on a Cu (111) substrate is un-der a nonuniform strain distribution. This nonuniformstrain suggests that there might be interesting PMFs inthe region of graphene surrounding the bubble. To an-alyze that we studied the PMF profile generated by theinflation of a graphene bubble constrained by a circularaperture with a radius of 4 nm on a Cu(111) substrate.The Cu-C interactions were modeled using a Morse po-tential with parameters D =0.1 eV, α =1.7 ˚A, r =2.2 ˚A,and a cutoff radius of 6 ˚A . Fig. 8 shows the PMF dis-tributions for differently shaped bubbles with deflectionof ∼ re-1 (a) (b) (c) (d) (e) (f) FIG. 8: (Color online) Top views of PMF patterns for graphene bubbles of different geometries on Cu (111) substrate. (a) circle(b) hexagon (c) pentagon (d) rectangle (aspect ratio 1:2) (e) square (f) triangle. All the bubble areas are ∼
50 nm , and sidelengths and pressures can be found in the main text. In all cases, the graphene lattice is oriented with the zig-zag directionalong the horizontal. The same color scale (in Tesla) is used in all panels. The edge of the substrate apertures used in the MDsimulations is outlined (gray line) for reference. flect a non-negligible graphene-Cu interaction, the PMFdistributions in Figs. 8(a-f) are much richer than inFigs. 6(a-f). That our simulation strategy involves press-ing graphene against the substrate certainly enhances theinteraction and promotes increased adhesion. This, inturn, adds a non-isotropic constraint for the longitudi-nal displacement and deformation of the graphene sheetwhich will affect the overall magnitude and spatial de-pendence of the PMF in the central region in such a waythat, for this case, the PMF magnitude is higher outsidethe inflated portion of graphene, rather than inside or inthe close vicinity of the boundary. This shows that thestrain and PMF patterns in graphene can be strongly in-fluenced by the chemical nature of the substrate and notjust its topography.To reveal the PMF that is induced by the substratealone, we show in Fig. 9 a side-by-side comparison ofthe PMFs that result when graphene is let to relaxon Au (111) and Cu (111), respectively. The plotteddata were obtained from energy minimization withoutpressure or aperture to show the intrinsic effect of thetwo substrates. Several interesting features emerge fromthese results, the first of which being the spontaneousdevelopment of a superlattice structure with a charac-teristic and well defined periodicity that is different in (a) (b) FIG. 9: (Color online) PMF distributions of graphene on per-fect (a) Cu (111) substrate and (b) Au (111) substrate withoutapertures nor gas pressure. The superlattice structure arisesnaturally from the need of the system to release strain buildupbecause of the mismatch in the lattice parameters of grapheneand the underlying substrate. The PMF scale is in units ofTesla. the two substrates. This Moir´e pattern in the PMF isthe result of a corresponding pattern in the strain fieldthroughout the graphene sheet, which is caused by theneed of the system to release strain buildup due to the2mismatch in the lattice parameters of graphene and thesubstrate. A second important aspect is the consider-able magnitude of the PMFs that can locally reach a fewhundreds of Tesla just by letting graphene reach the min-imum energy configuration in contact with the flat metalsubstrate. Another detail clearly illustrated by these twoexamples is the sensitivity to the details of the substrateinteraction: the substrate-induced PMF on Cu can bemany times larger than that on Au, and the Moir´e periodis also different. These super-periodicities are expected toperturb the intrinsic electronic structure of flat graphenewhose electrons now feel the influence of this additionalperiodic potential. That leads, for example, to the ap-pearance of band gaps at the edges of the folded Bril-louin zone. Such effects are currently a topic of interestin the context of transport and spectroscopic propertiesof graphene deposited on boron nitride, where this typeof epitaxial strain is conjectured to play a crucial role indetermining the metallic or insulator character .Since Fig. 9 reveals a strong graphene-substrate in-teraction, it is not surprising that the PMF patternsin Fig. 8 are still strongly dominated by the substrate-induced PMF. Unlike the cases discussed in Fig. 4, a sig-nificant structure remains in the PMF distribution out-side the hole region due to the tendency of the latticeto relax towards the characteristic Moir´e periodicity ofFig. 9(a) when in contact with a flat portion of substrate.In contrast, Au (111) has a larger lattice spacing and gen-erates considerably less epitaxial strain in the graphenefilm, implying comparatively weaker PMFs. It is thennatural that in the presence of the nanobubbles the ge-ometry of the aperture dominates the final PMF distribu-tion over the entire system when pressed against Au (111)(Fig. 6), whereas for Cu (111) the epitaxial contributionis the one that dominates (Fig. 8).
VI. BENDING EFFECTS
The large deflection-to-linear dimension ratio in theinflated graphene bubbles analyzed so far calls for ananalysis of the relative importance of the contribution tothe PMF from bending in comparison with that from thelocal stretching of the distance between carbon atoms.When full account of stretching and bending is takenby replacing the hopping (4) in the definition of thevector potential A given in Eq. (7) the resulting PMFcan have considerably higher magnitudes, as was alreadyseen in Fig. 2(f). To isolate the effect of bending aloneone can split the full hopping (4) in two contributions, t ij = t ( xy ) ij + t ( c ) ij , where the “in plane” stretching term issimply − t ( xy ) ij = V ppπ ( d ) . (10)Since the gauge field A is a linear function of the hopping(7), it can be likewise split into the respective stretch-ing and bending contributions so that A = A ( xy ) + A ( c ) . FIG. 10: (Color online) Density plot of the bending contribu-tion to the pseudomagnetic field, B ( c ) , for a circular graphenebubble with radius of 4 nm and a deflection of ∼ max(Bbc)/max(Bxy) R ( n m )
FIG. 11: (Color online) Ratio of the maximum PMF inducedby bending and stretching ( B c /B xy ) for circular graphenebubbles as a function of the graphene radius R , accordingto Hencky’s solution. When the PMF associated with A ( c ) is thus calculated forthe circular bubble of Fig. 2 we obtain the result shownin Fig. 10(a). As was already seen when comparing thedifferent PMF curves in Fig. 3, the effect of the curva-ture at the edges is quite remarkable and overwhelminglydominant in that region.More importantly, this fact could have been underap-preciated if the stretching and bending contributions hadbeen extracted only on the basis of an analytical solutionof the elastic problem such as Hencky’s model. To be def-inite in this regard let us consider the magnitude of thecontribution to the PMF that comes from bending in thecontinuum limit. If a gradient expansion of the full hop-ping (4) is performed, the vector potential (7) can beexpressed in terms of quadratic combinations of the sec-ond derivatives of the deflection h ( x, y ) . For example,3the term V ppπ ( d ) n i · n j in (4) leads to A ( c ) x = − a V ppπ qv F (cid:20)(cid:16) ∂ h∂y (cid:17) − (cid:16) ∂ h∂x (cid:17) (cid:21) , (11a) A ( c ) y = − a V ppπ qv F (cid:20) ∂ h∂x∂y (cid:16) ∂ h∂y + ∂ h∂x (cid:17)(cid:21) . (11b)This particular contribution was previously discussedby Kim and Castro Neto and, since all the bendingterms have the same scaling ∼ a h /R , where h and R are the characteristic height and radius, respectively,consideration of this one alone suffices for our purpose ofestablishing the magnitude of the bending terms in com-parison with the stretching one. Replacing the deflection h ( x, y ) provided by Hencky’s solution in Eqs. (11) leadsto the result shown in Fig. 10(b); it is clear that themaximum B c so obtained at the edges is much smallerthan the one derived from the atomistic simulation withthe full hopping. It is not surprising that the PMF com-ing from bending at the level of Hencky’s model is sosmall. A simple scaling analysis of the vector potentialsin the continuum limit shows that, from Eq. (8), A xy scales with strain as A xy ∼ (cid:15) and strain itself scales withdeflection as (cid:15) ∼ ( h/R ) for a characteristic linear dimen-sion R of the bubble. On the other hand, from (11) A c scales like A c ∼ ( ah ) /R . Therefore, the ratio B c /B xy will scale as ∼ ( a/R ) . Since the bubble under analy-sis has a/R ≈ .
04 the bending contribution is indeedexpected to be much smaller than the stretching one.We can even be more quantitative and extract the max-imum values of B c and B xy from Hencky’s solution andcompare their relative magnitudes as a function of circleradius, as shown in Fig. 11. Hencky’s solution predictsthat only when the radius of the circular bubble decreasesbelow about 1 nm does the contribution of the curvature-induced pseudomagnetic field become of the same orderas that due to in-plane stretching. This situation is equiv-alent to the need to account for the curvature and orbitalre-hybridization when describing the electronic structureof carbon nanotubes with diameters below length scalesof this same magnitude at the tight-binding level ;the neglect of these effects in the nanotube case leads toincorrect estimation of the band gaps and even of theirmetallic or insulating character.The problem with these considerations is that they failto anticipate the large effect at the edges, particularlythe scaling analysis which tells us only about the relativemagnitude of bending vs stretching in the central region.But, because we are inflating graphene under very highpressures in order to achieve deflections of the order of1 nm, a sharp bend results at the edge of the substrateaperture through which graphene can bulge outwards; itis this curvature effect that dominates the PMF plot inFig. 10, not the overall curvature of the bubble on thelarge scale. Hencky’s solution cannot capture this sinceit is built assuming zero radial bending moment at theedge . Moreover, since this happens within a distanceof the order of the lattice constant itself, the details of FIG. 12: (Color online) Density plot of the bending contribu-tion to the pseudomagnetic field, B ( c ) , for a graphene bubbledeflected to ∼ the displacements at the atomistic level including non-linearity and softening at large strains and curvaturesbecome crucial. This further highlights the importance ofaccurate atomistic descriptions of the deformation fieldsin small structures such as the sub-5 nm graphene bub-bles we have considered in this paper, and which havebeen shown experimentally to lead to significant pseudo-magnetic fields ; at this level models based on contin-uum elasticity theory can become increasingly limited foraccurate quantitative predictions and should be appliedwith caution.Finally, when realistic substrate conditions are con-sidered, one can see that the slippage effects contributevery differently for the PMFs arising from stretching andfrom bending. A general feature of the PMF distribu-tion obtained with realistic Au and Cu substrates is itssmaller overall magnitude in comparison with the artifi-cially clamped nanobubbles. This is easy to understandbecause the ability to slide in contact with the substrateallows graphene to stretch not only in the bubble region,but essentially everywhere, thereby reducing the strainconcentration around the edge of the aperture; and withsmaller strain gradients one gets smaller PMFs. Thebending effects, on the other hand, are not expected to bemuch affected by the sliding, especially when comparingnanobubbles with the same amount of vertical deflection,because the sharpness of the bend at the edge of theaperture is constrained mostly by the geometry alone.Direct inspection of the contribution to the PMF aris-ing from curvature in the Au and Cu substrates directlyconfirms this intuitive expectation, as shown in Fig. 12.Just as in the clamped case where graphene is pinnedto the substrate and cannot slide, the PMF associatedwith bending is seen to dominate the field distribution,with magnitudes similar to the registered in Fig. 10, andmuch larger than the PMF in the center of the bubble orthe substrate region (cf. Figs. 6 and 8). This not onlyshows how crucial the PMF associated with bending can4be in certain approaches to generate graphene nanobub-bles, but also that it is an effect largely insensitive to thedetails of the substrate. VII. DISCUSSION AND CONCLUSIONS
We have evaluated the strain-induced pseudomagneticfields in pressure-inflated graphene nanobubbles of differ-ent geometries and on different substrates whose config-urations under pressure were obtained by classical MDsimulations. The geometry of the nanobubbles is estab-lished by an aperture of prescribed shape in the substrateagainst which a graphene monolayer is pressed under gaspressure. Our results provide new insights into the na-ture of pseudomagnetic fields determined by the interplayof the bubble shape and the degree of interaction withthe underlying substrate. On a technical level, if bend-ing is (or can be) neglected, we have established that anapproximate displacement-based approach is adequate toobtain the strain tensor and accurate values of the pseu-domagnetic fields from MD simulations when comparedwith a direct tight-binding approach where the modifiedhoppings are considered explicitly.By comparing nanobubbles inflated in three differentsubstrate scenarios — namely, an arguably artificial, sim-ply clamped graphene sheet with no substrate couplingand more realistic conditions where the full interactionwith Au (111) and Cu (111) substrates is included fromthe outset in the MD simulations — we demonstratedthat the graphene-substrate interaction is an essentialaspect in determining the overall distribution and mag-nitude of strain and the PMFs both inside and outsidethe aperture region. For example, sections IV and Vdemonstrate that graphene can adhere substantially tothe substrate in atomically flat regions leading to sizablePMFs stemming only from epitaxial strain, even in theabsence of any pressure or substrate patterning. Thisadhesion varies from substrate to substrate and, in thepresence of an aperture or other substrate patterning,perturbs the final strain distribution of the nanobubblewhen compared with a simply clamped edge. On a morequantitative level, in the cases analyzed here where theaspect ratio of the bubbles is close to 1, the magnitudeof the PMFs associated with epitaxial strain alone caneasily be of the same magnitude as the PMF generatedwithin the bubble region. For Cu this is clear in Figs. 8and 9, and implies that the presence of the aperture isnot the main factor determining the field distribution.To better appreciate this aspect, we can inspect theaveraged cross-section of the PMF provided in Fig. 15whose the details are given in Appendix B. The key mes-sage conveyed by the data there is that under more real-istic conditions describing the graphene-substrate inter-action, and for the range of parameters explored here,the PMFs are no longer concentrated in and around theaperture. The section shown reveals that the PMF canbe considerably higher in regions well outside the aper- ture than inside or near the edge. This arises because, onthe one hand, the graphene-substrate interaction alone isable to generate considerable local strain gradients thatbeget PMFs as large as the ones that appear by forc-ing the inflation of graphene through the aperture [cf.Fig. 9]. On the other hand, the fact that graphene canslip into the aperture when simulated on the Au andCu substrates softens the strain gradients in its vicin-ity in comparison with the artificially clamped scenario.Slippage under these realistic substrate conditions pre-vents strain from concentrating solely within the aper-ture region which, instead, spreads to distances signifi-cantly away from the aperture. This is a sensible out-come on account of the very large stretching modulus ofgraphene that tends to penalize stretching as much aspossible. We illustrate this behavior in Fig. 13(a) thatcompares the magnitude of the radial displacement of cir-cular nanobubbles. Whereas in the artificially clampednanobubble graphene remains undisturbed (by design)outside the aperture, it is clear that in either the Au orCu substrates the carbon atoms pertaining to the regioninitially outside the aperture are radially pulled every-where towards it under pressure, as one intuitively ex-pects. One consequence of this is the softening of straingradients in the bubble region: slippage naturally tendsto diminish the PMFs generated within and around theaperture. The other is that, obviously, the deflection atthe center is increased, as shown in Fig. 13(b).The joint effect of these two factors (slippage and ad-hesion to the substrate) is that forcing graphene intoa nanobubble profile at the center of the system is nolonger effective in concentrating the strain gradients and,consequently, the PMF is no longer more prominentthere. One immediate implication of this is the fact thatwhether or not it is feasible to locally tailor the PMF dis-tribution on very small (nanometer) scales depends notonly on the elastic response of graphene or its loadingand geometric constraints, but also on the nature of thesubstrate involved.It is also clear from the above that there might existcertain substrates in which the epitaxial strain can be sig-nificant enough to, by itself, lead to visible modificationsof the electronic structure of graphene, and even lead tomodified transport characteristics . Incidentally ourpressure-based approach facilitates and promotes a uni-form adhesion because graphene is compressed againstthe substrate. It would be interesting to experimentallystudy graphene on top of such substrates inside pressurechambers, and assess the degree of control that can beachieved over the Moire patterns and the modificationsof the electronic and transport characteristics.Another important factor to consider in estimating themagnitude and profile of the PMF generated under agiven set of force distributions and geometric constraintsis whether those conditions lead to strong local curva-ture of the graphene lattice. We analyzed this issue hereby separately considering the contributions from bondbending and from stretching to the PMF in the repre-5 (a) (b)
FIG. 13: (Color online) The radial and vertical displacementcomponents of graphene along the diametral section y = 0of three circular nanobubbles of radius 4 nm, correspondingto the three substrate conditions considered in this report.Panel (a) shows the component u x of the graphene in-planedisplacement which corresponds, in essence, to the radial dis-placement because of the circular symmetry. Panel (b) showsthe vertical deflection. All cases were inflated under the samepressure of ≈
19 kBar. sentative case of circular nanobubbles. Our results es-tablish that, even though the overall, large-scale curva-ture of the graphene sheet leads to significant correctionsto the pseudomagnetic field only in ultra-small bubbleswith diameter smaller than 2 nm, sharp bends arisingfrom direct clamping or from being pressed against anedge in the substrate aperture result in much strongerPMFs locally. At the qualitative level this is naturallyexpected and certainly not surprising. What is surprisingsignificant is that the bending contribution can be manytimes larger than its stretching counterpart, leading toa PMF distribution dominated by large values near theedges of the substrate apertures. Moreover, since this isa local geometric effect, it does not depend on the bubblesize but only on the local curvature around sharp bends,and should remain in considerably larger systems. Thisindicates that curvature of the graphene sheet should cer-tainly not be ignored in many situations involving out-of-plane deflection, even though the scaling analysis basedon the overall profile could point otherwise. Finally, we underline once more that the strategy togenerate graphene nanobulges through gas pressure waschosen here to minimize other external forces and influ-ences on the deflection and slippage of graphene whilebeing able to produce deflections and aspect ratios equalto those reported in recent experiments that explore thelocal electronic properties of these structures. But theconclusions and implications discussed above certainlycarry to various other means of achieving such or similarnanobubbles and have, therefore, a wide reach and wideimport beyond graphene pressurized through apertures.
Note added . Recently, we became aware of a recentproposal to connect structure and electronic properties oftwo-dimensional crystals based on concepts from discretegeometry that allows yet another efficient alternative toobtain the strain and PMF at discrete lattice points with-out the need, for example, to perform numerical deriva-tives upon the displacement fields or vector potentialsextracted from the MD data . Acknowledgments
ZQ acknowledges support from a Boston UniversityDean’s Catalyst Award and the Mechanical EngineeringDepartment at Boston University. ZQ and HSP acknowl-edge support of U.S. National Science Foundation grantCMMI-1036460. The authors thank Prof. Teng Li forilluminating discussions on the stress method. VMP andAHCN are supported by the NRF CRP grant “Novel 2Dmaterials with tailored properties: beyond graphene” (R-144-000-295-281). This work was supported in part bythe U.S. National Science Foundation under grant No.PHYS-1066293 and the hospitality of the Aspen Centerfor Physics.
Appendix A: Angular averaging of the PMF
Fig. 14 below shows the radial dependence of the av-eraged PMF amplitude close to the edge of the circularaperture for the clamped circular case discussed in sec-tion III (Fig. 2). In Fig. 3 we plot the amplitude of thePMF at the edge of the circular aperture in the substratefor various inflation pressures with clamped graphene. InFig. 5(b) we show the average amplitude of the PMF atdifferent distances from the center.In all these cases, the data shown reflect the PMF am-plitude averaged over the azimuthal direction. To extractthe average PMF at a given radius, the 2D distributionof the field is divided into a sequence of radial and az-imuthal bins (annular sectors). For each radial annulusthere are 20 bins, each with a 18 degree width. Thewidth of the radial annulus is chosen such that at least10 atoms lie in each bin (this is why there are fewer datapoints near the center of the bubble). The average andstandard deviations of the PMF in each bin correspondto the value and error bar of that bin. For example, each6 (cid:72) Å (cid:76) P M F (cid:144) S i n (cid:72) Θ (cid:76)(cid:72) T (cid:76) TB, completeDispTB, in (cid:45) planeHencky
FIG. 14: (Color online) Angular-averaged amplitude of thePMF for the same cases presented in Fig. 2 in the form ofdensity plots. The horizontal axis represents the distancefrom the center of a pressurized circular graphene nanobub-ble with clamped boundary conditions. The data containedhere is the same shown in Fig. 3(a), except that here the (or-ange) data corresponding to the PMF obtained from the fullhopping perturbation [Eq. (4)] is included for comparison aswell. The bending effects are clearly dominant around theedge/clamping region. Away from the edge, and inside, thethree numerical curves follow Hencky’s model.
20 40 60 80 - - - - y H Å L P M F H T L ClampedGoldCopper
FIG. 15: (Color online) The PMF of graphene pressurizedthrough equally sized circular apertures along a vertical sec-tion extending from the center of the aperture to the bottomof the simulation cell. The vertical line at ∼
40 ˚A marks theradius of the apertures. See the text in Appendix B for moredetails. point in Fig. 5(c) corresponds to this average PMF for agiven bin. Afterwards, for each radial annulus the datais fit to the expected sin(3 θ ) dependence. The amplitudeof the best fit is plotted as a point [e.g., as in Fig. 14]and the fitting error provides the error bar. Appendix B: Sectional plot of PMFs for circularapertures
To better illustrate the magnitude of the PMF in thevicinity of apertures simulated under realistic substrateconditions, we present here a sectional view of the fieldfor the representative cases of the circular apertures sim-ulated in the artificially clamped, Au, and Cu scenariosexplored in the main text. Fig. 15 shows the PMF ofgraphene pressurized through circular apertures of thesame size sampled along a vertical section extending fromthe center of the aperture to the bottom of the simula-tion cell. The sections are taken from the correspondingdata shown in Figs. 4(a), 6(a), and 8(a) by sampling thePMF along a vertical direction and performing averageswithin square bins of 25 ˚A . The averaging is done to ac-count for local fluctuations in the PMF and the standarddeviation in each bin is used to draw the error bars. Thetraces in Fig. 15 are analogous to the ones in Figs. 14 or3(a), with the exception that there is no angular aver-aging here because the substrate interaction breaks therotational symmetry [cf., for example, the region outsidethe aperture in Fig. 8]; consequently the error bars arehigher here than in the artificially clamped cases. Appendix C: Comparison of PMFs fromdisplacement and full TB approaches
As described in the main text, the displacement ap-proach to obtain the pseudomagnetic fields through-out the graphene sheet consists in directly employingEq. (8), where the components of the strain tensor areextracted numerically from the MD-relaxed atomic posi-tions. Apart from contributions beyond linear order instrain, this should be equivalent to computing the vectorpotential A ( R ) directly from the definitions (7) and (4),but neglecting the bending effects in the hopping. Thisamounts to considering − t ij = V ppπ ( d ).For completeness, and to show that the two approacheslead to the same results in practice, we present in Fig. 16the PMF distribution computed by the displacement ap-proach for the same systems analyzed in Fig. 4. Theagreement is very satisfactory and shows that the dis-placement and tight-binding methods are equivalent ifcurvature can be neglected. Appendix D: Comparison of Displacement andStress Approaches
The final (inflated bubble) configuration gives us thebasic ingredients needed to calculate the strain, i.e. , thedeformed atomic positions. Here we present further de-tails on the displacement and stress approaches we inves-tigated for calculating the strain. In the end, the stressapproach revealed itself inadequate to accurately capturethe local strain in the graphene lattice.7
FIG. 16: (Color online) PMF distribution for the same systems analyzed in Fig. 4, but here the field is computed by thedisplacement approach discussed in the main text. The PMF scale is in units of Tesla. The edge of the substrate aperturesused in the MD simulations is outlined (gray line) for reference.
1. Displacement Approach
We begin with the continuum definition for strain ,which is written as (cid:15) ij = 12 (cid:18) ∂u i ∂X j + ∂u j ∂X i (cid:19) + 12 ∂u k ∂X i ∂u k ∂X j , (D1)where (cid:15) ij are the components of the strain, u is the dis-placement, and X denotes the position of a point in thereference configuration. To compute the displacementfield that is needed to evaluate the strain in Eq. D1,we first exploit the geometry of the graphene latticeby meshing it using tetrahedral finite elements , whereeach finite element is comprised of four atoms. To re-move spurious rigid body translation and rotation modes,we choose the deformed position of the atom of interest(atom 0) to be the reference position, i.e. , r = R .By subtracting the original position of each neigh-boring atom from its deformed position, we obtain thedisplacement vectors of the three nearest neighbors: u = ( u x , u y , u z ), u = ( u x , u y , u z ), u =( u x , u y , u z ).We use the linear interpolation property of the four-node tetrahedral element to denote the displacement field U ( x, y, z ) = ( U x , U y , U z ) inside the tetrahedral elementas: U x = a x + a y + a z , U y = a x + a y + a z , U z = a x + a y + a z , where a to a are unknown constants for each tetrahedral element. Inserting the positions( r = ( x , y , z ), r = ( x , y , z ), r = ( x , y , z ))and the corresponding displacements ( u , u , u ) ofthe three neighboring atoms, we can solve a to a interms of r , r , r and u , u and u , thus obtainingall coefficients of U ( x, y, z ). If we rearrange U ( x, y, z ) toexpress it in terms of u , u and u , we obtain thefollowing equation: U x U y U z = N N N N N N
00 0 N N N u x u y u z u x u y u z u x u y u z , (D2)where N i = N i ( x, y, z ) , i = 1 , , U = N · u N , (D3)where u N = [ u , u , u ] T is the displacement field ofthe three neighbor atoms.After we obtain the displacement field U , the straincan be derived by differentiating Eq. (D3) following the8continuum strain as defined in Eq. (D1) to give (cid:15) = T · u N , (D4)where T = ∂ N ∂ x is constant inside each tetrahedral ele-ment. Once the strains for each atom are determined, thevector gauge field A is straightforward to compute. How-ever, to get the pseudomagnetic field, B = ∂ x A y − ∂ y A x ,another derivative is needed, calculated in a similar fash-ion as the strain is calculated from the displacement field.Thus, the displacement approach involves two numericalderivatives, but no approximation is made about mate-rial properties.
2. Stress Approach
In MD simulations, the atomic virial stress can be ex-tracted on a per-atom basis. In the present work, thevirial stress as calculated from LAMMPS was obtainedfor the final (inflated) graphene bubble configuration.These stresses were then related to the strain via a linearconstitutive relationship, as was done recently by Klimovet al. . In the current work, we utilized a plane stressmodel for graphene, where the in-plane strains are writ-ten as (cid:15) xx = E ( σ xx − µσ yy ), (cid:15) yy = E ( σ yy − µσ xx ), (cid:15) xy = σ xy G . The material properties of graphene are cho-sen as E = 1 TPa , G = 0.47 TPa and µ = 0.165 ,where E is the Young’s modulus, G the shear modulus,and µ Poisson’s ratio. It is important to note that, be-cause a linear stress-strain relationship is assumed, theresulting strain is generally underestimated , particularlyat large deformations due to the well-known nonlinearstress-strain response of graphene .Both potential and kinetic parts were taken intoaccount for virial stress calculation. We note thatthe virial stress calculated in LAMMPS is in units of“Pressure · Volume”, and thus we used the standardvalue of 3.42 ˚A as the effective thickness of single layergraphene to calculate the stress. A plane stress consti-tutive model was utilized to calculate the strain via (cid:15) xx (cid:15) yy (cid:15) xy = E − µE − µE E
00 0 G · σ xx σ yy σ xy , (D5)where the constitutive parameters are given in the maintext of the manuscript.After the strain is obtained, the same method as inthe displacement approach was used to calculate the vec-tor gauge field A and the pseudomagnetic field B . Thestress approach avoids one numerical differentiation buta constitutive approximation is involved, i.e. , that thestress-strain response for graphene is always linear.
3. Benchmark Examples
We compare the displacement and stress approachesvia two simple benchmark examples, those of uniaxial stretching and simple shear. For the uniaxial stretchingcase, (cid:15) xx ≈
10 % strain was applied along the x-direction.The loading was done by applying a ramp displacementthat went from zero in the middle of simulation box to amaximum value at the +x and -x edges of the graphenemonolayers.For the simple shear case, (cid:15) xy ≈ code withthe AIREBO potential . The result for the uniaxialstretching is shown in Figs. 17 and Fig. 18, while thesimple shear is shown in Figs. 19 and Fig. 20. The supe-rior performance of the displacement approach is seen inboth cases. Specifically, because a linear stress-strain re-lationship is assumed in the stress approach as shown inEq. (D5), the resulting strain is generally underestimated ,particularly at large deformations due to the well-knownnonlinear stress-strain response of graphene .Once the strain distribution is determined from theMD simulations the PMF, B , can be directly evaluatedfrom the definitions above. However, if the strain tensoris calculated within the deformation approach, a secondnumerical derivative is needed to get B , which is likelyto introduce a certain degree of error. Nevertheless wefound the errors to be of acceptable magnitude.Compared with the displacement approach, the stressapproach avoids one numerical differentiation, but a con-stitutive approximation is involved. To compare the ac-curacy of the displacement and stress approaches, we cal-culated the PMF distribution in a circular bubble (forwhich an analytic solution is available and detailed anal-ysis was recently performed ) by obtaining the strainvia three different methods, as illustrated in Fig. 2: ananalytic continuum mechanics model, i.e. , the Henckysolution (b), the MD-based displacement approach (c),and the MD-based stress approach (d). In the MD sim-ulations we used 100 snapshots over 5 ps during thermalequilibrium to determine the average final position andstress for the inflated bubbles. For all three models, theradius of the circular hole was 3 nm, while the final de-flection was about 1 nm.As Fig. 2 demonstrates, the PMFs generated from theMD-based displacement approach are in good agreementwith those that follow from Hencky’s analytic solution,and also with previously reported values for a circularbubble . In contrast, the stress approach fails to yieldreasonable results for this loading situation, even at thequalitative level.9 x(angstrom) (cid:161) xx distribution y ( ang s t r o m ) FIG. 17: (Color online) (cid:15) xx distribution by displacement ap-proach for uniaxial stretching case with 10 % strain. x(angstrom) (cid:161) xx distribution y ( ang s t r o m ) FIG. 18: (Color online) (cid:15) xx distribution by stress approachfor uniaxial stretching case with 10 % strain. x(angstrom) (cid:161) xy distribution y ( ang s t r o m ) − − − − − FIG. 19: (Color online) (cid:15) xy distribution by displacement ap-proach for simple shear case with 1 % strain. x(angstrom) (cid:161) xy distribution y ( ang s t r o m ) − − − − − FIG. 20: (Color online) (cid:15) xy distribution by stress approachfor simple shear case with 1 % strain. § Corresponding author: [email protected] ∗∗ Corresponding author: [email protected] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, andA. A. Firsov, Nature , 197 (2005). A. K. Geim and K. S. Novoselov, Nature Materials , 183(2007). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Reviews of Modern Physics , 109 (2009). M. Y. Han, B. Ozyilmaz, Y. Zhang, and P. Kim, PhysicalReview Letters , 206805 (2007). J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken,M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, et al.,Science , 213 (2010). F. Guinea, M. I. Katsnelson, and A. K. Geim, NaturePhysics , 30 (2010). Z. Qi, D. A. Bahamon, V. M. Pereira, H. S. Park, D. K.Campbell, and A. H. Castro Neto, Nano Letters , 2692(2013). H. Tomori, A. Kanda, H. Goto, Y. Ootuka, K. Tsukagoshi,S. Moriyama, E. Watanabe, and D. Tsuya, Applied PhysicsExpress , 3 (2011). F. Guinea, A. K. Geim, M. I. Katsnelson, and K. S.Novoselov, Physical Review B , 035408 (2010). V. M. Pereira and A. H. Castro Neto, Physical ReviewLetters , 046801 (2009). F. Guinea and T. Low, Philosophical Transactions of theRoyal Society a-Mathematical Physical and EngineeringSciences , 5391 (2010). F. Guinea, B. Horovitz, and P. Le Doussal, Physical Re-view B , 205421 (2008). N. Abedpour, R. Asgari, and F. Guinea, Physical ReviewB , 115437 (2011). K.-J. Kim, Y. M. Blanter, and K.-H. Ahn, Physical ReviewB , 081401 (2011). N. C. Yeh, M. L. Teague, S. Yeom, B. L. Standley, R. T. P.Wu, D. A. Boyd, and M. W. Bockrath, Surface Science , 1649 (2011). H. T. Yang, Journal of Physics-Condensed Matter (2011). A. L. Kitt, V. M. Pereira, A. K. Swan, and B. B. Goldberg,Physical Review B , 115432 (2012). A. L. Kitt, V. M. Pereira, A. K. Swan, and B. B. Goldberg,Physical Review B , 159909(E) (2013). K. Yue, W. Gao, R. Huang, and K. M. Liechti, Journal of Applied Physics , 083512 (2012). V. M. Pereira, A. H. Castro Neto, H. Y. Liang, and L. Ma-hadevan, Physical Review Letters , 156603 (2010). D. A. Abanin and D. A. Pesin, Physical Review Letters , 066802 (2012). Z. F. Wang, Y. Zhang, and F. Liu, Physical Review B ,041403 (2011). V. M. Pereira, A. H. Castro Neto, and N. M. R. Peres,Physical Review B , 045401 (2009). Z. H. Ni, T. Yu, Y. H. Lu, Y. Y. Wang, Y. P. Feng, andZ. X. Shen, ACS Nano , 483 (2009). M. Farjam and H. Rafii-Tabar, Physical Review B ,167401 (2009). S.-M. Choi, S.-H. Jhi, and Y.-W. Son, Physical Review B , 081407(R)+ (2010). C. L. Kane and E. J. Mele, Physical Review Letters ,1932 (1997). H. Suzuura and T. Ando, Physical Review B , 235412(2002). V. M. Pereira, R. M. Ribeiro, N. M. R. Peres, and A. H.Castro Neto, Eur. Phys. Lett. , 67001 (2010). L. Gong, I. A. Kinloch, R. J. Young, I. Riaz, R. Jalil, andK. S. Novoselov, Adv. Mater. , 2694 (2010). N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui,A. Zettl, F. Guinea, A. H. Castro Neto, and M. F. Crom-mie, Science , 544 (2010). J. Lu, A. H. Castro Neto, and K. P. Loh, Nature commu-nications , 823 (2012). M. Neek-Amal and F. M. Peeters, Physical Review B ,195446 (2012). M. Neek-Amal, L. Covaci, and F. M. Peeters, PhysicalReview B , 041405 (2012). W. Humphrey, A. Dalke, and K. Schulten, Journal ofMolecular Graphics , 33 (1996). Lammps, http://lammps.sandia.gov (2012). S. Plimpton, Journal of Computational Physics , 1(1995). S. J. Stuart, A. B. Tutein, and J. A. Harrison, Journal ofChemical Physics , 6472 (2000). H. Zhao and N. R. Aluru, Journal of Applied Physics ,064321 (2010). Z. Qi, F. Zhao, X. Zhou, Z. Sun, H. S. Park, and H. Wu,Nanotechnology , 265702 (2010). M. Neek-Amal and F. M. Peeters, Physical Review B ,195445 (2012). W. G. Hoover, Physical Review A , 1695 (1985). R. R. Nair, H. A. Wu, P. N. Jayaram, I. V. Grigorieva, andA. K. Geim, Science , 442 (2012). J. S. Bunch, S. S. Verbridge, J. S. Alden, A. M. van derZande, J. M. Parpia, H. G. Craighead, and P. L. McEuen,Nano Letters , 2458 (2008). S. P. Koenig, N. G. Boddeti, M. L. Dunn, and J. S. Bunch,Nature Nanotechnology , 543 (2011). N. N. Klimov, S. Jung, S. Z. Zhu, T. Li, C. A. Wright, S. D.Solares, D. B. Newell, N. B. Zhitenev, and J. A. Stroscio,Science , 1557 (2012). T. J. R. Hughes,
The Finite Element Method: LinearStatic and Dynamic Finite Element Analysis (Prentice-Hall, 1987). A. Isacsson, L. M. Jonsson, J. M. Kinaret, and M. Jonson,Physical Review B , 035423 (2008). A. Hansson and S. Stafstr¨om, Physical Review B ,075406 (2003). F. de Juan, J. L. Ma˜nes, and M. A. H. Vozmediano, Phys- ical Review B , 165131 (2013). J. V. Sloan, A. Sanjuan, Z. Wang, C. Horvath, andS. Barraza-Lopez, Physical Review B , 155436 (2013). M. Oliva-Leyva and G. G. Naumis, Physical Review B ,085430 (2013). P. Wang, W. Gao, Z. Cao, K. M. Liechti, and R. Huang,Journal of Applied Mechanics , 040905 (2013). W. B. Fichter, Tech. Rep. 3658, NASA Langley ResearchCenter (1997). Y. Wei, B. Wang, J. Wu, R. Yang, and M. L. Dunn, NanoLetters , 26 (2013). C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science ,385 (2008). S. Jun, T. Tashi, and H. S. Park, Journal of Nanomaterials , 380286 (2011). A. De Martino, L. Dell’Anna, and R. Egger, Physical Re-view Letters , 066802 (2007). M. R. Masir, P. Vasilopoulos, and F. M. Peeters, New J.Phys. , 095009 (2009). C. H. Y. X. Lim, A. Sorkin, Q. Bao, A. Li, K. Zhang,M. Nesladek, and K. P. Loh, Nature Communications ,1556 (2012). A. Kitt, Z. Qi, S. Remi, H. S. Park, A. K. Swan, andB. Goldberg, Nano Letters , 2605 (2013). S. Piana and A. Bilic, Journal of Physical Chemistry B , 23467 (2006). R. E. Tuzun, D. W. Noid, B. G. Sumpter, and R. C.Merkle, Nanotechnology , 241 (1996). A. Rytknen, S. Valkealahti, and M. Manninen, Journal ofChemical Physics , 5826 (1998). R. He, L. Zhao, N. Petrone, K. S. Kim, M. Roth, J. Hone,P. Kim, A. Pasupathy, and A. Pinczuk, Nano Letters ,2408 (2012). K. Maekawa and A. Itoh, Wear , 115 (1995), ISSN0043-1648. C. Woods, L. Britnell, A. Eckmann, R. Ma, J. Lu, H. Guo,X. Lin, G. Yu, Y. Cao, R. Gorbachev, et al., NaturePhysics , 451 (2014). P. San-Jose, A. Guti´errez-Rubio, M. Sturla, and F. Guinea,Phys. Rev. B , 075428 (2014). M. Neek-Amal and F. M. Peeters, Applied Physics Letters , 041909 (2014). J. Jung, A. DaSilva, S. Adam, and A. H. MacDonald(2014), arXiv:1403.0496. E. A. Kim and A. H. Castro Neto, Europhysics Letters (2008). X. Blase, L. X. Benedict, E. L. Shirley, and S. G. Louie,Physical Review Letters , 1878 (1994). S. Barraza-Lopez, A. A. Pacheco Sanjuan, Z. Wang, andM. Vanevi´c, Solid State Communications , 70 (2013). A. A. Pacheco Sanjuan, Z. Wang, H. P. Imani, M. Vanevi´c,and S. Barraza-Lopez, Phys. Rev. B (2014). J. A. Zimmerman, D. J. Bammann, and H. J. Gao, Inter-national Journal of Solids and Structures , 238 (2009). A. P. Thompson, S. J. Plimpton, and W. Mattson, Journalof Chemical Physics (2009). Y. Huang, J. Wu, and K. C. Hwang, Physical Review B , 245413 (2006). K. Min and N. R. Aluru, Applied Physics Letters ,013113 (2011). O. L. Blakslee, D. G. Proctor, E. J. Seldin, G. B. Spence,and T. Weng, Journal of Applied Physics41