Relating elasticity and graphene folding conformation

Abstract

Variational calculus is employed to determine the folding behaviour of a single graphene sheet. Both the elastic and van der Waals energies are taken into account, and from these considerations the shape of the curve is determined. By prescribing that the separation distance between the folded graphene in the parallel region is 3.32 Å, an arbitrary constant arising by integrating the Euler–Lagrange equation is determined, and the full parametric representations for the folding conformation are derived. Using typical values of the bending rigidity in the range of 0.800–1.60 eV, the shortest stable folded graphene sheets are required to be at least 6.5–10 nm in length.

---

1 Introduction

Carbon atoms on a graphene sheet are connected as a planar hexagonal array of sp²-bonds. This two-dimensional system is one of the most promising materials for nanoelectromechanical systems,^1–3 and may also be used in biological applications^4–6 because its bending stiffness is comparable with that of lipid bilayers.⁷ Two graphene layers connected by a continuous folded curved graphene edge exhibit a number of interesting electronic,^8–12 magnetic,¹³ thermal¹⁴ and mechanical^7,15 properties. Moreover, it has been shown that self-folding graphene is stable and driven by van der Waals interactions.^16–18 Here, by taking into account both elastic and van der Waals energies, the calculus of variations is exploited to study the self-folding property of a graphene sheet or graphitic nano ribbon (GNR).

Cranford et al.¹⁹ employed a coarse-grained model where pseudo-atoms were used to represent groups of atoms to study the folding behaviour for single- and multi-layers of graphene, and found that this behaviour depends on the bending rigidity and the surface energy of the sheets. Both analytical and numerical methods have been carried out by Rainis et al.¹³ to determine the deformation of folded graphene. Further, Lopez-Bezanilla et al.¹² used a density functional theory calculation to determine the geometric structure of a single- and multi-layer of folded graphene. This theoretical finding was also confirmed by their experimental images taken using high-resolution transmission electron microscopy, and these authors suggested that the number of graphene layers is a factor to the folding profile. Similarly, Meng et al.¹⁸ undertook molecular dynamics simulations and beam theory to investigate the self-folding of a single-layer graphene, and found that the chirality of the folded edge does not influence the folding behaviour. Graphene aerogels have high- and super-compressible properties due to the elasticity of the individual sheets and they have potential applications for vibration and shock damping. Further, graphene folds are the elementary constituents of graphene aerogels responsible for the observed unusual compressibility of aerogels.²⁰

In formulating a first approximation of graphene folding following,²¹ we assume that graphitic materials deform as in perfect elasticity, and rather like the elastica.²¹ As a consequence, we propose an elastic bending energy as being proportional to the square of the curvature and we employ variational calculus to determine the folded conformation which minimises the bending energy while simultaneously maximising the van der Waals interaction for the folded graphene sheet, noting that both the curvature effect arising from the elastic deformation and the inter-spacing between two graphene layers arising from the van der Waals energy are accommodated. Although this analytical model is simple, involving only two parameters which are the van der Waals interaction strength and the bending rigidity, it successfully describes the folding conformation of graphene sheets. We comment that the same variational approach to minimise the elastic energy has been adopted by the present authors in a number of joining problems involving carbon nanostructures. The present approach is developed from ref. 21–24 and we refer the reader to the review paper ref. 25 which includes further references.

In the following section, the proposed mathematical model is formulated and the variational formulations are described in Sections 3 and 4. Further, a full parametric solution involving the bending rigidity is given in Section 5, and finally some concluding remarks are made in Section 6.

2 Conformation model

Here we propose to model the curve formed by a single sheet of graphene folding over onto itself. In this configuration we suggest that the two dominant forces are those originating from the elastic folding deformation of the graphene sheet and the van der Waals interactions between the sheet. With reference to Fig. 1, we assume a rectangular graphene sheet for which one edge folds uniformly over onto the sheet and that the consequent folding geometry both symmetric with respect to the x-axis and two-dimensional in the sense that there are no variations in the perpendicular z-direction. Furthermore, assuming a reflective symmetry in the folded geometry, we need only consider the top half of the problem. The shown curve may be divided into three parts, which are denoted by C_n where n = {1, 2, 3}. C₁ is the curve from (0, 0) to (x₀, y₀) where the line curvature of the fold is negative, C₂ is the curve from (x₀, y₀) to (x₁, δ) where the corresponding line curvature is positive and C₃ is the straight line from (x₁, δ) to (x₂, δ) where the line curvature is zero. We will denote the concatenation of these three regions as C = C₁ + C₂ + C₃.

Fig. 1 Geometry for folded graphene with total half length L.

We assume that the elastic energy E_e may be calculated by integrating the square of the curvature κ over the length of the total curve multiplied by a scaling constant γ, which is the bending rigidity of graphene. Therefore, we may write

where κ denotes the line curvature of the curve y = y(x) and is given by

We now assume that the region dominated by van der Waals interactions comprises only the C₃ portion of the solution curve and we model the van der Waals energy E_v using the formulation

where is a positive constant giving the van der Waals interaction energy per area length and u(x − x₁) is the Heaviside unit step function. We now consider the total energy normalised by the magnitude of the bending rigidity γ which we denote by E, therefore
where we have introduced the constant , the ratio of van der Waals interaction strength to the bending rigidity, both of which are material characteristics of graphene and therefore prescribed.

3 Variational formulation

We now employ calculus of variations to determine the shape of the curve C₁ + C₂, treating the quantity E as a functional of a conformation profile y = y(x) and the particular solution that we seek is the one for which E takes the minimum value. We also impose an isoperimetric constraint on the total arclength of C, namely
which we incorporate into the variational principle using a Lagrange multiplier λ. Therefore, the functional under consideration here is given by
Noting that both the curvature effect and the van der Waals contribution are taken into account.

Now considering the arclength constraint we comment that for C₃ then y′ = 0 and ds = dx, and therefore

and hence

We note that L is a curved half length of the graphene sheet. Moreover, we comment that u(x − x₁) = 0 for C₁ + C₂ and κ = 0 for C₃ then

We further comment here that F depends linearly on the Lagrange multiplier λ and it takes no part in determining the shape of y = y(x) which minimises the value of F. Therefore, we are justified in choosing any value for λ and for convenience we choose λ = α. Hence, the functional we wish to minimise is given by

subject to the endpoint conditions y(0) = 0, and y′(0) → ∞, and at x = x₁ we have a natural boundary condition on x which for a second order variational problem requires thatwhere the subscripts denote partial differentiation and f is the integrand of the functional F{y} given by

Also on the endpoint x = x₁, we have that y(x₁) = δ and y′(x₁) = 0.

4 Solution to the variational problem

Since the function contains no explicit dependence on y, we may integrate the Euler–Lagrange equation once to obtainwhere β is an arbitrary constant of integration. Furthermore, the function is also independent of x, and therefore we have
where H is a constant. However, we know from the natural boundary condition (3) that H must equal zero at x = x₁ and since H is a constant then H = 0 throughout the domain, and therefore we may deduce
We also note that the term in parentheses is precisely the first integral of the Euler–Lagrange equation given in (5) and it is a constant, so that the solution must satisfy

f − βy′ − y′′f_y′′ = 0.

Now on substituting f(y′, y′′) given in (4) into this equation, we obtain
and after some rearrangement, we may deducewhere again κ is the line curvature defined by (1).

We now seek a parametric solution using the substitution y′ = tan θ to obtain

and hence parametrically we obtain the first order equations
For convenience, we now make the substitution θ = 2ϕ − π/2, and after some rearrangement we have

The differential equation for x(ϕ) given in (7)₁ can be integrated immediately giving

where c₁ is an arbitrary constant.

The differential equation for y(ϕ) given in (7)₂ may be rearranged as follows

where the terms involving ϕ can be integrated as elliptic integrals of the first and second kinds, therefore we havewhere c₂ is an arbitrary constant of integration, F and E denote the incomplete elliptic integrals of the first and second kinds, respectively, and the elliptic modulus k = [2β/(α + β)]^1/2.

5 Full parametric solution

Critical to the parametric solution is the point (x₀, y₀), where the curvature changes sign and it is also the boundary point between the two curves C₁ and C₂. This point corresponds to the parameter value ϕ = ϕ₀, where ϕ₀ = sin⁻¹(1/k).

On the curve C₁, the solution varies over the range ϕ = [0, ϕ₀] starting at the origin and having a trajectory that is upward and bending to the right half of the plane. To satisfy these requirements we require

where the starting point (0, 0) has been used to determine that c₁ = 2(α + β)^1/2/β and c₂ = 0. This means the critical point (x₀, y₀) is given by

For C₂ we change the signs adopted from the general solutions (8) and (9), and we require that the solution for C₁ and C₂ are continuous at (x₀, y₀), which leads to the solution given by

where ϕ ∈ [π/4, ϕ₀], and again F and E are the incomplete elliptic integrals of the first and second kinds, respectively.

The parameter value α corresponds to the material characteristics where it is defined as the ratio of the van der Waals interaction , to the bending rigidity γ. For the graphene, it was estimated by Spanu et al.²⁶ that ≈ 0.0214 eV Å⁻². Following the work by Wei et al.,⁷ the possible values of γ range from 0.800 to 1.60 eV. Therefore, we have the corresponding values of α in the range of 0.0268 to 0.0134 Å⁻². In fact, the shape function or the folding conformation is a function of γ, and it can depend on the number of defects and the number of dopings. Indirectly, the folding conformation could therefore indicate the presence of defects and dopings.

To determine the remaining unknown β, we use the endpoint , where δ is a half of the separation distance between the folded graphene in the parallel region, and we adopt the value δ ≈ 1.66 Å.²⁷ Solving this equation numerically, the value of β can be determined. A plot of five possible folded graphenes is shown in Fig. 2, and the numerical values are presented in Table 1.

Fig. 2 Fold region for various values of bending rigidity γ.

Table 1 Numerical values for graphene folding using van der Waals interaction to be 0.0214 eV Å⁻² (ref. 26)

γ (eV)	x0 (Å)	y0 (Å)	Lmin (Å)
0.800	7.88	2.50	13.7
1.00	9.19	2.74	16.2
1.20	10.3	2.95	18.3
1.40	11.4	3.15	20.2
1.60	12.3	3.32	22.0

---

By integrating ds over both of the curves C₁ + C₂, we may determine the arclength of the extremal curve and this provides a lower bound to L for this solution to be physically meaningful, and we may deduce

where dx/dϕ and dy/dϕ are given in (7). Finally, we havewhere L_min is a lower bound to the half length L.

The total energy for the folded graphene is

and on changing the integration variable using y′ = tan θ, we may deduce

Again, we use θ = 2ϕ − π/2 so that E_tot can be evaluated in terms of elliptic integrals of the second kind and the final expression becomes

We note that E(0, k) = 0 and (x₂ − x₁) = L − L_min. Fig. 3 shows the relation between the total energy E_tot of the folded graphene and the total half length L, and the system is stable provided that the total energy is negative. For a bending rigidity γ = 1.40 eV, we require that the graphene sheet have a half length of at least 45 Å for the folding behaviour to occur and this is in agreement with Meng et al.¹⁸

Fig. 3 Relation between total energy E_tot and total half length L where L for L ≥ L_min.

Furthermore, the folding conformations obtained by this model are also adopted to compare with the folding images taken by Lopez-Bezanilla et al.¹² and they are as shown in Fig. 4. We observe that the best fit profiles shown in Fig. 4 (top) and (bottom) are obtained by using different bending rigidity values which are γ = 1.00 eV (navy colour) and γ = 1.60 eV (purple colour), respectively. Note that the model presented here is derived only for a single layer of graphene, and the different values for γ are due to the increased bending rigidity for multi-layers. Moreover, the five folding conformations in each figure are obtained independently but they are plotted on the same coordinates for the purpose of comparison. Overall, the folding geometry of graphene as determined here by the variational calculus provides a simple model but correctly describes the folding profile observed in experiments.

Fig. 4 Folding conformations obtained by present model superimposed upon the images taken by high-resolution transmission electron microscopy¹² where the proposed model perfectly matches with (top) γ = 1.00 eV represented by the navy line and (bottom) γ = 1.60 eV represented by the purple line (adapted with permission from Lopez-Bezanilla et al.¹² Copyright 2012 American Chemical Society).

6 Summary

In this paper, we use a calculus of variations approach to determine the shape of folded graphene sheets. The calculus of variations is utilised to minimise the elastic energy arising from the curvature squared while maximising the van der Waals energy in the parallel region. In addition, the fold is assumed to have a translational symmetry along the fold, so that the problem may be reduced to a two dimensional problem with reflective symmetry across the fold. Here, the van der Waals interaction is assumed to be 0.0214 eV Å⁻² and with a separation distance in the parallel region fixed to be 1.66 Å. However, the bending rigidity of graphene is selected in the range from 0.800 to 1.60 eV. An analytical parametric solution is obtained and numerical values for the critical points where the curvature changes the sign and where the folded graphene assumes a parallel orientation are determined. Furthermore, the model predicts the critical length of a graphene sheet needed to adopt an energetically stable folded conformation that depends on the actual value of the bending rigidity adopted for the graphene sheet.

References

1. J. S. Bunch, A. M. van der Zande, S. S. Verbridge, I. W. Frank, D. M. Tanenbaum, J. M. Parpia, H. G. Craighead and P. L. McEuen, Electromechanical resonators from graphene sheets, Science, 2007, 315, 490–493.
2. K.-T. Lam, C. Lee and G. Liang, Bilayer graphene nanoribbon nanoelectromechanical system device: a computational study, Appl. Phys. Lett., 2009, 95, 143107.
3. M. Poetschke, C. G. Rocha, L. E. F. Foa Torres, S. Roche and G. Cuniberti, Modeling graphene-based nanoelectromechanical devices, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 193404.
4. X. Sun, Z. Liu, K. Welsher, J. T. Robinson, A. Goodwin, S. Zaric and H. Dai, Nano-graphene oxide for cellular imaging and drug delivery, Nano Res., 2008, 1, 203–212.
5. H. Jiang, Chemical preparation of graphene-based nanomaterials and their applications in chemical and biological sensors, Small, 2011, 7, 2413–2427.
6. Y. Wang, Z. Li, J. Wang, J. Li and Y. Lin, Graphene and graphene oxide: biofunctionalization and applications in biotechnology, Trends Biotechnol., 2011, 29, 205–212.
7. Y. Wei, B. Wang, J. Wu, R. Yang and M. L. Dunn, Bending rigidity and Gaussian bending stiffness of single-layered graphene, Nano Lett., 2013, 13, 26–30.
8. K. Nakada, M. Fujita, G. Dresselhaus and M. S. Dresselhaus, Edge state in graphene ribbons: Nanometer size effect and edge shape dependence, Phys. Rev. B: Condens. Matter, 1996, 54, 17954.
9. T. Enokia, Y. Kobayashia and K.-I. Fukuia, Electronic structures of graphene edges and nanographene, Int. Rev. Phys. Chem., 2007, 26, 609–645.
10. K. A. Ritter and J. W. Lyding, The influence of edge structure on the electronic properties of graphene quantum dots and nanoribbons, Nat. Mater., 2009, 8, 235–242.
11. J. Feng, L. Qi, J. Y. Huang and J. Li, Geometric and electronic structure of graphene bilayer edges, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 165407.
12. A. Lopez-Bezanilla, J. Campos-Delgado, B. G. Sumpter, D. L. Baptista, T. Hayashi, Y. A. Kim, H. Muramatsu, M. Endo, C. A. Achete, M. Terrones and V. Meunier, Geometric and electronic structure of closed graphene edges, J. Phys. Chem. Lett., 2012, 3, 2097–2102.
13. D. Rainis, F. Taddei, M. Polini, G. Leon, F. Guinea and V. I. Falko, Gauge fields and interferometry in folded graphene, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 165403.
14. K. V. Zakharchenko, M. I. Katsnelson and A. Fasolino, Finite temperature lattice properties of graphene beyond the quasiharmonic approximation, Phys. Rev. Lett., 2009, 102, 046808.
15. C. Lee, X. Wei, J. W. Kysar and J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene, Science, 2008, 321, 385–388.
16. H. C. Schniepp, K. N. Kudin, J.-L. Li, R. K. Prudhomme, R. Car, D. A. Saville and I. A. Aksay, Bending properties of single functionalized graphene sheets probed by atomic force microscopy, ACS Nano, 2008, 2, 2577–2584.
17. J. Zhang, J. Xiao, X. Meng, C. Monroe, Y. Huang and J.-M. Zuo, Free folding of suspended graphene sheets by random mechanical stimulation, Phys. Rev. Lett., 2010, 104, 166805.
18. X. Meng, M. Li, Z. Kang, X. Zhang and J. Xiao, Mechanics of self-folding of single-layer graphene, J. Phys. D: Appl. Phys., 2013, 46, 055308.
19. S. Cranford, D. Sen and M. J. Buehler, Meso-origami: Folding multilayer graphene sheets, Appl. Phys. Lett., 2009, 95, 123121.
20. H. Hu, Z. Zhao, W. Wan, Y. Gogotsi and J. Qiu, Ultralight and highly compressible graphene aerogels, Adv. Mater., 2013, 25, 2219–2223.
21. J. Zang, A. Treibergs, Y. Han and F. Liu, Geometric constant defining shape transitions of carbon nanotubes under pressure, Phys. Rev. Lett., 2004, 92, 105501.
22. C. Li and T.-W. Chou, A structural mechanics approach for the analysis of carbon nanotubes, Int. J. Solids Struct., 2003, 40, 2487–2499.
23. Y. Jin and F. G. Yuan, Simulation of elastic properties of single-walled carbon nanotubes, Compos. Sci. Technol., 2003, 63, 1507–1515.
24. T. Natsuki, K. Tantrakarn and M. Endo, Effects of carbon nanotube structures on mechanical properties, Appl. Phys. A: Mater. Sci. Process., 2004, 79, 117–124.
25. D. Baowan, B. J. Cox and J. M. Hill, Determination of join region between carbon nanostructures using variational calculus, ANZIAM J., 2013, 54, 221–247.
26. L. Spanu, S. Sorella and G. Galli, Nature and strength of interlayer binding in graphite, Phys. Rev. Lett., 2009, 103, 196401.
27. L. A. Girifalco, M. Hodak and R. S. Lee, Carbon nanotubes, buckyballs, ropes, and a universal graphitic potential, Phys. Rev. B: Condens. Matter, 2000, 62, 13104–13110.

---