Limits of stability in supported graphene nanoribbons subject to bending

Graphene nanoribbons are prone to in-plane bending even when supported on flat substrates. However, the amount of bending that ribbons can stably withstand remains poorly known. Here, by using molecular dynamics simulations, we study the stability limits of 0.5-1.9 nm wide armchair and zigzag graphene nanoribbons subject to bending. We observe that the limits for maximum stable curvatures are below ~10 deg/nm, in case the bending is externally forced and the limit is caused by buckling instability. Furthermore, it turns out that the limits for maximum stable curvatures are also below ~10 deg/nm, in case the bending is not forced and the limit arises only from the corrugated potential energy landscape due to the substrate. Both of the stability limits lower rapidly when ribbons widen. These results agree with recent experiments and can be understood by means of transparent elasticity models.

Today graphene nanoribbons can be fabricated at atomic precision but only in the presence of a stabilizing substrate [1]. The substrate stabilizes flimsy ribbons and suppresses their tendency to twist, fold, and ripple [2][3][4][5][6]. However, even substrates cannot fully prevent all deformations, most of which induce mechanical strains that alter ribbons' electronic properties [7][8][9]. Actually, such strain engineering of electronic properties is gaining popularity, whereby detailed knowledge of mechanical stability limits is becoming increasingly valuable [10].
Mechanical strain can be created, for example, by lattice mismatch, by impurities, by lattice defects, and by the fabrication process itself [11]. Compressive strain, in particular, is often limited by buckling instability. For uniaxial compression buckling has been observed in experiments at 0.5% strain and in simulations at 0.8% strain [12,13]. In graphene nanoribbons, however, the most pertinent deformation is not uniaxial compression but bending. Yet, the mechanical stability limits of supported ribbons subject to bending remain unexplored. In this paper, therefore, we aimed to address two fundamental questions: How much can a graphene nanoribbon of a given width bend on a given substrate until it buckles? And, to what extent can it remain bent due to the corrugation potential energy of the substrate alone without external forcing? As it will turn out, both of these questions could be answered by transparent modeling.
Our simulations were closely related to the recent experiments of van der Lit et al. in Ref. [14] (Fig. 1). There an atomically precise seven-armchair graphene nanoribbon was bent at low temperatures on a Au(111) surface by an atomic force microscope (AFM) tip. Under forced bending and above certain maximum curvature the ribbon was observed to buckle off the substrate [ Fig. 1(a)]. Furthermore, the ribbon was observed to withstand certain maximum curvature, presumably due to the lateral energy corrugations arising solely from the substrate interactions [ Fig. 1(c)].
To investigate the buckling instability in more detail, we simulated ribbons subject to forced bending [ Fig. 1(b)]. We * pekka.koskinen@iki.fi simulated hydrogen-passivated N -armchair (N = 5, 7, 9, 11, and 13) and N -zigzag (N = 4, 6, 8, and 10) graphene nanoribbons of widths w ≈ 0.5-1.9 nm and lengths given by 1/10 aspect ratio. The C-C, C-H, and H-H interactions were modeled by the reactive empirical bond-order (REBO) potential [15]. The ribbons were initially relaxed on a model Au substrate, which assumed an interaction with the ribbon described by a z-dependent potential with 20-meV/Å 2 adhesion, 3.4-Å equilibrium distance, and a functional form suggested FIG. 1. Seven-armchair graphene nanoribbons subject to bending. (a) In experiments bending was controlled by the tip of an atomic force microscope (AFM), whose movements are denoted by arrows. Buckling is seen as the bright kink. (b) In simulations ribbons were bent by fixing their front ends and by turning their tail ends. The rightmost geometry shows the buckled geometry. (c) Maximum curvature without external forcing. After manipulation the ribbon remained bent by the substrate corrugations alone. Scale bar, 10 nm. (d) In the simulations one end of the ribbon (green tail) was pinned to (set in registry with) the substrate while the other end was turned to the maximum stable curvature beyond which the entire ribbon started sliding. The experimental figures in panels (a) and (c) are reproduced from Ref. [14] by Creative Commons Attribution licence; image ordering has been changed. by the Lennard-Jones (LJ) 12-6 potential ( Fig. 2) [16][17][18]. This substrate model ignores lateral energy corrugation, but it is expected to be a good approximation because graphene nanoribbons that are out of registry with respect to the Au(111) substrate have been shown not to experience any lateral forces and thus to exhibit superlubricity [19,20].
The supported ribbons were simulated by the LAMMPS code, using a 1-fs time step and a Langevin thermostat at 10-K temperature and 5-ps damping time [22]. First the ribbons were thermalized on the model substrate. Then they were gradually bent by fixing one end and slowly (quasistatically) turning the other end while simultaneously allowing its free movement on the plane [ Fig. 1(b)]. At a later instant the turning direction was reversed, and the simulation terminated with straight ribbons.
At the initial stages of the simulations the bending was smooth, and the ribbons remained adhered to the substrate. Here we quantify the amount of bending both by the in-plane curvature κ = 1/R, where R is the radius of curvature, and by the dimensionless curvature = κw/2, which also equals the absolute amount of strain at the ribbon edges. Using straightforward continuum elasticity theory, the elastic energy during this initial stage is where w is the ribbon width, l is the ribbon length, k = 19 eV/Å 2 is graphene's in-plane modulus, and τ is the stress at the passivated armchair (τ ac = −1.5 eV/Å) or zigzag edges (τ zz = −0.2 eV/Å) as given by the REBO potential [23].
During this initial stage we observed weak ripples at the inner edges of the ac ribbons. Ripples were notable up and down displacements of alternating armchair units and they were observable along the entire ribbon. They have been observed also in straight ribbons where they have been attributed to chemically induced edge stress; here the edge stress was created mostly by the bent geometry itself [24,25]. When curvature increased, the rippling amplitude increased, but the wavelength remained fixed. These ripples were observed only for the ac ribbons as the zz ribbons remained almost completely flat prior to buckling.
When the increasing curvature reached a critical limit, the in-plane stress finally became unbearable, and the ribbon suddenly buckled [ Fig. 1(c)]. Buckling allowed two parts of the ribbon to straighten, which released in-plane elastic energy, although at the expense of lost adhesion and increased out-of-plane bending energy. Buckling occurred later for narrow ribbons than for wide ribbons. The events during the bending-straightening simulations are best gauged through the maximum height of the ribbon above the substrate [ Fig. 3(b)]. Initially the buckle was formed at b , but upon straightening it remained stable also for curvatures < b so that when the ribbon finally unbuckled at b , roughly half the buckling curvature, the result was a notable hysteresis. The buckling-unbuckling process was reversible; plastic deformations did not occur. These observations are in agreement with experiments that also showed the restoring of the initial geometry. In particular, for the N = 7 ac ribbon the buckling occurred in experiments at a curvature of 4 deg/nm, in reasonable agreement with the computational curvature of 6 deg/nm [14]. Note that it is justifiable to compare experiments only to the smaller curvature b because in macroscopic time scales random perturbations help drive the system toward buckled geometry already at smaller curvatures.
To understand the general width dependence in the buckling [ Fig. 3(c)], let us develop a model that accounts for the in-plane strain, out-of-plane bending, and substrate adhesion energies. In the model the ribbon is treated as two aligned narrow strands that represent the compressed and stretched halves of the ribbon. The aligned strands are next to the neutral line and separated by w eff = αw where the widthdependent parameter α( 1) is later fitted to account for the averaging. Upon buckling the outer strand remains flat, but the height profile of the inner strand acquires the form y(l) = A sin 2 (l/λπ) (0 l λ), where A is the buckling amplitude and l is the distance measured along the strand. This profile decreases the strand length by l = π 2 A 2 /(4λ) and thereby relieves the compressive strain energy at the inner edge by wk l/2 and the tensile strain energy at the outer edge by the same amount. This approach is similar to that in Ref. [26]. Adding this strain energy release to the loss in Lennard-Jones energy [ w/2[V LJ (y) − V LJ (u)]dl] and the out-of-plane bending energy associated with the height profile [ w 4 Dy 2 (l)dl], the energy difference between purely bent and buckled ribbon becomes Here vdw is the adhesion energy per unit area, σ is the interlayer distance, and D = 1.0 eV is graphene's bending modulus [26,27]. Buckling occurs when the first term becomes large enough due to the increasing curvature so that E( b ) = 0. The energy of the buckled geometry is further minimized by ∂ E/∂λ| λ=λ b = 0. Solving these equations yields λ b = 9Å and Fit to the simulations gives α i = 1/(β i w −1 + 1), where β ac = 7 and β zz = 5Å, which provide a good agreement with the simulated buckling curvatures [ Fig. 3(c)]. The fit is physically meaningful and obeys the consistency requirement α 1. The validity of the model is probably limited for ribbon widths below a few nanometers, although b = 2.3% is a reasonable limit for very wide ribbons too. Although our simulations included ribbons only with hydrogen-passivated zigzag and armchair edges, other edges with other passivations or edge reconstructions are possible [28][29][30][31]. Especially in free-standing graphene the edges may create sizable corrugations [25,32]. On substrates these corrugations diminish in magnitude but do not vanish completely [33]. However, here the edge stresses are small due to hydrogen passivation and the lateral stresses due to bending are so large that the effect of edge stress is fairly small. This is suggested already by the quantitatively similar buckling behavior in zigzag and armchair ribbons [ Fig. 3(c)].
For completeness, we repeated buckling simulations for armchair ribbons also at room temperature. As the main result, the effect of temperature was to reduce the hysteresis and initiate buckling at slightly smaller curvatures [ Fig. 3(c)]. On average, however, the buckling occurred around the same curvature as described by the model fitted at low temperatures.
In the next set of simulations, we investigated the limits of maximal curvature in armchair ribbons allowed by the substrate energy corrugation alone. In these simulations we chose to place the ribbons on a graphene substrate modeled by the Kolmogorov-Crespi (KC) registry-dependent interlayer potential [21]. This model substrate was obviously different from the Au(111) substrate in the experiments, but our choice was a necessary compromise for a feasible substrate model with a realistic energy corrugation. Namely, the frequently used Lennard-Jones potential typically yields an order of magnitude with too low-energy corrugation for sliding, and the proper registry-dependent potentials for graphene and Au(111) are missing [34]. Nevertheless, the ribbon adhesions for both Au and graphene substrates are similar, so KC potential was an attempt to combine a well-defined substrate model with a realistic corrugation energy landscape.
In these simulations one end of the ribbon was first appended by a tail of length L t that was pinned to the substrate by setting it in full registry [ Fig. 1(f)]. The other end was then gradually turned until the maximum stable curvature beyond which the pinning was released and the tail started sliding, causing straightening of the ribbon. The ribbon was considered stable at given L t and κ if it remained in place for 20 ps, although it was evident already within a few picoseconds whether the curvature was stable or not. The maximum curvature limits were then searched for each ribbon width with several tail lengths.
Simulations show that narrow ribbons withstand higher curvatures than wide ribbons and that maximum curvatures increase when the tail lengths increase (Fig. 4). It is notable that certain finite curvatures can be achieved even in the absence of any added tail (the inset of Fig. 4). This occurs because the ribbon's end also is close to registry and not yet subject to superlubric behavior. By geometry considerations we therefore approximate that the length L t = √ 2Ra + a 2 close to the end of the ribbon is still pinned to the substrate, where a is a length scale for the tolerance in a lateral displacement that is still considered to be in registry. Thus, the total length of the substrate-pinned ribbon at the end equals L pin = L t + L t . This assumption serves as a starting point for a model for the maximum curvature limit. In the model we consider the pinned part to be subject to a bending moment kw 2 /6 imposed by the unpinned part. This moment must not exceed a maximum value, lest the pinned part starts to slide. At the maximum the bending moment equals the maximum allowed moment, or where f is the maximum force per unit area during sliding, averaged over all sliding directions. The integration is over the pinned part of length L pin , and r is the distance to its center of mass. Fitting the force parameter f with a chosen tolerance a = 0.7Å to simulation data yields f = 0.7 meV/Å 3 . The maximal force per unit area for sliding in an armchair direction is f max = 2.3 meV/Å 3 , which confirms the physical interpretation of the fit (f ≈ 0.3f max ) [35]. Upon inserting these parameters into the model Eq. (4), the trends in maximum curvatures get reproduced surprisingly well (Fig. 4). The model underestimates the maximum curvatures as compared to simulations, which is however not surprising given the highly discrete nature of the short-tail limit (inset of Fig. 4). The model predicts pinning at roughly constant edge strain of ∼0.9%, but in simulations the allowed edge strain depends somewhat on ribbon width, changing as ribbons widen from ∼0.9% for N = 5 to ∼1.5% for N = 13. Such dependence may originate due to thermal fluctuations, which affect narrow and wide ribbons differently due to the different number of pinned atoms.
These simulations can be compared to the experimentally observed pinning in Ref. [14], although with caution. The energy corrugations for graphene ribbons on Au(111) and on graphene are probably different but likely of similar magnitude due to the similarity of the adhesion itself [16]. To this end, note that the model in Eq. (4) suggests that the substrate affects the trends only through the averaged parameter f . Thus, even though the symmetry in Au(111) differs from that in graphene, it is not unreasonable to expect that the results would correspond also to Au substrate, at least semiquantitatively. Such correspondence is further supported by the rough agreement between the experimental (2 deg/nm for the gold substrate) and the simulated (1.4 deg/nm for the model graphene substrate) maximum curvatures for a seven-armchair ribbon [14]. At any rate, the parameter f allows transferring the results to any other substrate, making the model highly versatile.
Whereas in buckling the effect of temperature was clearly small, in pinning its effect is more ambiguous. Although the energy corrugation per atom ∼ 9 meV corresponds only to the temperature of T ∼ 100 K, the pinning also still occurred at room temperature, at least within time scales accessible to the simulations (20 ps). The general tendency of an increased temperature was to modestly decrease the maximal pinning curvature, although the results became less clear. Whereas at low temperatures the possible unpinning of the tail was fast (∼1 to 2 ps) and clear cut, at high temperatures thermal fluctuations brought unambiguity by introducing more variations to the time scale of unpinning. Thus, reliable determination of structure stability would have required simulation times beyond reasonable limits as also indicated by recently observed sliding phenomena [36].
To conclude, these simulations and the associated models provide transparent understanding for the stability limits in supported graphene nanoribbons subject to bending. Narrow 5-, 7-, 9-, 11-, and 13-armchair ribbons require only minimal pinned parts to maintain curvatures around 1 deg/nm (radius of curvature R ≈ 60 nm). Although such curvatures are gentle, other studies have found them to cause predictable modifications in ribbons' electronic and optical properties. In particular, the simulations in Ref. [37] showed that the energy gap for N -armchair graphene nanoribbons changes according to the expression, where q = mod(N,3) (restricted to q = 0,1) is the ribbon family, γ = 1.7 describes bond anharmonicity that is relevant for bending-induced stretching, and δ = 12 eV is an electromechanical coupling constant related to gap changes during the stretching of straight ribbons. Combining Eq. (5) with Eq. (3), the buckling-limited maximum energy gap change becomes directly E max g (w) = 5.4 × (7Å × w −1 + 1) 2 meV.
For the ribbons studied here this amounts from 23-meV (N = 5) to 10-meV (N = 13) gap changes. For wider ribbons the maximum gap change shrinks. In the case of pinning the maximum curvature depends on the tail length L t , but it is always limited by Eq. (3), so with unconstrained bending, Eq. (6) gives the upper limit for gap changes. Buckling, however, can modify the electronic properties even more than bending. Simulations showed that narrow ribbons remained flat above 4-deg/nm curvatures (R ≈ 14 nm), but stability was strongly width dependent; ribbons wider than 1.5 nm remained flat only below 2 deg/nm (R 29 nm). The obtained stability limits thus provide guidelines to design experiments and to choose structures that would be stable enough for reliable device operation. Because the adhesion energies for most van der Waals bound physisorbed twodimensional materials are of similar magnitude, we expect the presented elastic models to have applicability for several other ribbon and substrate materials [17]. To this end, we propose that the stabilities of bent ribbons could even be used as a measurement technique to investigate the interaction between different nanoribbons and substrates.
We thank the Academy of Finland for funding (Projects No. 283103 and No. 251216) and the CSC-IT Center for Science in Finland for computer resources.