The Gaussian Stiffness of Graphene deduced from a Continuum Model based on Molecular Dynamics Potentials
Cesare Davini Antonino Favata Roberto Paroni

	















 The Gaussian Stiffness of  Graphene 
  deduced from a Continuum Model 
 based on Molecular Dynamics Potentials
    
Cesare Davini^1 Antonino Favata^2 Roberto Paroni^3

    December 30, 2023
================================================================================================================







	
	^1 Via Parenzo 17, 33100 Udine

		mailto:cesare.davini@uniud.itcesare.davini@uniud.it

		
		^2 Department of Structural and Geotechnical Engineering

		Sapienza University of Rome, Rome, Italy

		mailto:antonino.favata@uniroma1.itantonino.favata@uniroma1.it

		
		
		^3 Dipartimento di Architettura, Design e Urbanistica

		University of Sassari, Alghero (SS), Italy

		mailto:paroni@uniss.itparoni@uniss.it
	





myheadings
C. Davini, A. Favata, R. Paroni 
The Bending Behavior of  Graphene





§ ABSTRACT


We consider a discrete model of a graphene sheet with atomic interactions governed by a harmonic approximation of the 2nd-generation Brenner potential that depends on bond lengths, bond angles, and two types of dihedral angles. A continuum limit is then deduced that fully describes the bending behavior. In particular, we deduce for the first time an analytical expression of the Gaussian stiffness, a scarcely investigated parameter  ruling the rippling of graphene,  for which contradictory values have been proposed in the literature. We disclose the atomic-scale sources of both bending and Gaussian stiffnesses and provide for them quantitative evaluations.



Keywords: Graphene, Continuum Modeling, Gaussian stiffness.









§ INTRODUCTION



Graphene   has attracted increasing interest during the past few years,  and is nowadays used in a great variety of applications, taking advantage of its extraordinary mechanical, electrical and thermal conductivity properties. Nevertheless, its potentialities, and those of graphene-based materials, are far from being fully explored and exploited, and many studies are carried out by the scientific community in order to develop new technological applications <cit.>.

The understanding of the bending behavior of graphene is  of paramount importance in several technological applications. It is exploited, for example, to predict the performance of graphene nano-electro-mechanical devices and ripple formation <cit.>, and it is proposed to be the key point to produce efficient hydrogen-storage devices <cit.>. 

A very recent review in Materials Today  by  Deng &  Berry <cit.> gives an overview on the hot problem of wrinkling, rippling and crumpling of graphene, highlighting formation mechanism and applications. Indeed, these corrugations  can modify its electronic structure, create polarized carrier puddles, induce pseudo-magnetic field in bilayers and alter surface properties. Although a great effort has been done on the experimental side, predictive models are still wanted. They are of crucial importance when these phenomena need to be controlled and designed.

In particular, since the bending stiffness and the Gaussian  stiffness —that is the reluctance to form non-null Gaussian curvatures— are the two crucial parameters  governing the rippling of  graphene, it is necessary to accurately determine them for both the design and the manipulation of graphene morphology. Although  several evaluations of the bending stiffness have been proposed in the literature, the Gaussian stiffness has not been object of an intensive study. Indeed, as pointed out in a very recent review on mechanical properties of graphene <cit.>, only two conflicting evaluations have been proposed. In <cit.>, periodic boundary conditions have been used within a quantum-mechanical framework, and the value -0.7 eV has been found. While  in <cit.> the estimate of -1.52 eV has been obtained by combining the configurational energy of membranes determined by Helfrich Hamiltonian with energies of fullerenes and single wall carbon nanotubes calculated by Density Functional Theory (DFT). At a discrete level there are  two main  difficulties in the evaluation of the Gaussian stiffness: on the one hand, controlling a discrete double curvature surfaces is problematic, and on the other hand, a suitable notion of Gaussian curvature at the discrete level should be introduced. Instead, when well established continuum models are adopted, such as plate theory,  one has the problem of determining the equivalent stiffnesses, letting alone the conceptual crux of giving a meaning to  the notion of thickness (see <cit.>, <cit.> and references therein).

In this paper we deduce  a continuum 2-dimensional model of a graphene sheet inferred from Molecular Dynamics (MD). In particular, looking at the 2nd-generation reactive empirical bond-order (REBO) potential <cit.>, we give a nano-scale description of the atomic interactions and then we deduce the continuum limit, avoiding the problem of postulating an “equivalent thickness” and circumventing artificial procedures to  identify the material parameters that describe the mechanical response of a  plate within the classical theory.

Our analysis of the atomic-scale interaction relies on the discrete mechanical model proposed in <cit.> and exploited in  <cit.>, whose results are also based on the  2nd-generation Brenner potential. This potential is largely used in MD simulations for carbon allotropes; for a detailed description of its general form  and that adopted in our theory we  refer the reader to  Appendix B of <cit.>. Here, we recall the key ingredients needed:

	
  * the kinematic variables associated with the interatomic bonds involve first, second and third nearest neighbors of any given atom. In particular, the kinematical variables we consider are bond lengths, bond angles, and dihedral angles; from <cit.> it results that these latter are of two kinds, that we here term C and Z, as carefully described in Sec. <ref>.
	
  * graphene suffers an angular self-stress,  and the self-energy associated with the self-stress (sometimes called cohesive energy in the literature) is quantitatively relevant;
	
  * the energetic contribution of dihedral interaction is very relevant in bending.


For the first time, we propose a continuum model able to predict both the bending  and the Gaussian stiffnesses. The analytical formula we obtain for the former predicts exactly the same value as that computed with MD simulations of the last generation. The value of the Gaussian stiffness we obtain is in very good agreement with DFT computations proposed in <cit.>. 

For the modeling of graphene many different approaches at different scales can be found in the literature, ranging from first principle calculations <cit.>, atomistic calculations <cit.> and continuum mechanics <cit.>. Furthermore, mixed atomistic formulations with finite elements have been reported for graphene <cit.>.


The paper is organized as follows. In Sec. <ref>, we describe the kinematics and the energetic of the  graphene sheet at the nano-scale. In Sec. <ref>, we deduce the  strain measures  for the change of edge lengths, wedge angles and dihedral angles, approximated to the lowest order that makes the energy quadratic in the displacement. In Sec. <ref>, the total energy is split in its in-plane and out-of-plane contributions, and focus is set on the latter, having the first already been considered in <cit.>. In Sec. <ref>, we deduce a continuous energy that approximates the discrete energy and in Sec. <ref> the limit energy is rearranged in a more amenable form, able to put in evidence the equivalence with plate theory.  In Sec. <ref>, quantitative results for the continuum material parameters are deduced, by means of the 2nd-generation Brenner potential, and  compared with the literature.  Appendix <ref>, containing some computations ancillary to Sec. <ref>, completes the paper.





§ DESCRIPTION OF KINEMATICS AND ENERGETICS OF THE GRAPHENE SHEET
 

At the nano-scale a graphene sheet is a discrete set of carbon atoms that, in the
absence of external forces, sit at the vertices of a
periodic array of hexagonal cells. More specifically,
atoms occupy the nodes of the 2–lattice, see Figure <ref>,  generated by two simple Bravais lattices


    [ L_1(ℓ) = {𝐱∈ℝ^2:
    𝐱 = n^1ℓ_1 + n^2ℓ_2     (n^1,
    n^2) ∈ℤ^2 },;                                       L_2(ℓ) = ℓ𝐩 + L_1(ℓ), ]

simply shifted with respect to one another. In (<ref>), ℓ denotes
the lattice size (the reference interatomic distance), while
ℓ_α and ℓ𝐩 respectively are the lattice
	vectors and the shift vector, whose Cartesian components
are given by

    _1 =(√(3), 0),   _2 = (√(3)/2,
    3/2 )     𝐩 =
    (√(3)/2, 1/2).


The sides of the hexagonal cells in Figure <ref> stand for the bonds
between pairs of next nearest neighbor atoms and are represented by the
vectors

    _α = _α -    (α = 1, 2)     _3 = -  .

For convenience we also set 
    _3=_2-_1.



As reference configuration we take the set of points ^ℓ∈ L^1(ℓ)∪ L^2(ℓ) contained in a bounded open set Ω⊂ℝ^2.

Graphene mechanics is ruled by the interactions between
the carbon atoms given by some suitable potential. According to the 2nd-generation Brenner potential <cit.>, as detailed in <cit.>, in order to account properly for the mechanical behavior of a bended graphene sheet it is necessary to consider three types of energetic contributions, respectively coming from: binary interactions between next nearest atoms  (edge bonds), three-bodies interactions between consecutive pairs of next nearest atoms  (wedge bonds) and four-bodies interactions between three consecutive pairs of next nearest atoms (dihedral bonds).  There are two types of relevant dihedral bonds: the Z-dihedra
in which the edges connecting the four atoms form a z-shape and the C-dihedra in which the edges form a c-shape (see Fig. <ref>).





We consider a harmonic approximation of the stored energy and assume that it is given by the sum of the following terms:


    [                      𝒰_ℓ^l= 1/2 ∑_ℰ k^l   (l  -
    l^)^2 ,;                      𝒰_ℓ^ϑ = 1/2 ∑_𝒲 k^ϑ  (ϑ  -  ϑ^)^2,; 𝒰_ℓ^Θ = 1/2 ∑_𝒵 k^  (  -  Θ^)^2+1/2 ∑_𝒞 k^  (  -  Θ^)^2 ]

𝒰_ℓ^l, 𝒰_ℓ^ϑ and 𝒰_ℓ^Θ are the energies of the edge bonds, the wedge bonds and the dihedral bonds, respectively; l denotes the distance between nearest neighbor atoms, ϑ the angle between pairs of edges having a lattice point in common and  and  the Z- and C-dihedral angles between two consecutive wedges, to be  defined later (see Fig. <ref>); l^ is the edge length at ease, 
ϑ^ the angle  at ease between consecutive edges and 
Θ^ the dihedral angle at ease. 




The sums extend to
all edges, ℰ, all wedges, 𝒲, all Z-dihedra, 𝒵, and all C-dihedra, , contained in the set Ω. The
bond constants k^l, k^ϑ, k^,  and k^  will be deduced by making use of the 2nd-generation Brenner potential.

The graphene sheet does not have a configuration at ease (i.e. stress-free). Indeed, in <cit.> it has been shown that 

    Θ^=0,    l^ =ℓ    ϑ^=2/3π + δϑ_0,

where δϑ_0 0.
We set δΘ:=Θ, l= ℓ + δ l and ϑ=2/3π + δϑ and write (<ref>) as

    [             𝒰_ℓ^l = 1/2 ∑_ℰ k^l   (δ l)^2 ,;       𝒰_ℓ^ϑ = 1/2 ∑_𝒲 k^ϑ  (δϑ  -  δϑ_0)^2,; 𝒰_ℓ^Θ = 1/2 ∑_𝒵 k^  (δ)^2+1/2 ∑_𝒞 k^  (δ)^2 ]

In particular, up to a constant, the wedge energy takes the form

    𝒰_ℓ^ϑ = τ_0∑_𝒲δϑ +  1/2 ∑_𝒲 k^ϑ  (δϑ)^2,

with 

    τ_0 := -k^ϑ δϑ_0
 
the angle self-stress.
The dihedral bonds play an important role because they contribute to the stored energy by about 50%, see <cit.>, the rest is due to the angle self-stress τ_0 associated to the wedge bonds. 

The energy decomposition (<ref>) is based on the choice of the set of kinematical variables {l, ϑ, , }. This choice is the most natural, if one considers the 2nd-generation Brenner potential, where all those variables  appear in explicit manner.  A harmonic approximation  in each of those parameters is of course unique. 	

In the next section we shall make explicit the change of length δ l, the change of wedge angle δϑ, and the changes of the Z- and C-dihedral angles δ and δ.
In Section <ref>, with the notation introduced in the next section, we shall write the energies
(<ref>) more explicitly.






§ APPROXIMATED STRAIN MEASURES

In this section we calculate the strain measures associated to a change of configuration described by a displacement field : (L_1(ℓ)∪ L_2(ℓ))∩Ω→ℝ^3, approximated to the lowest order that makes the energy quadratic in .



 §.§ Change of the edge lengths

With δ l_i(^ℓ) we denote the change in length of the edge parallel to _i and starting from the lattice point ^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω. Thus,

    δ l_i(^ℓ)   =|(^ℓ + ℓ_i + (^ℓ+ℓ_i)) - (^ℓ +(^ℓ))|-ℓ
       =|ℓ_i + ((^ℓ+ℓ_i) - (^ℓ))|-ℓ,

and up to terms o(| |) can be rewritten as

    δ l_i (^ℓ)= 1/ℓ ((^ℓ+ℓ_i) - (^ℓ))·𝐩_i      i=1,2,3.

In particular, the first order changes are determined by the in-plane components of  only.



 §.§ Change of the wedge angles

For each fixed node ^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω we denote by ϑ_i(^ℓ) the angle of the wedge delimited by the edges _i+1 and _i+2; that is, the wedge angle opposite to the i-th edge (see Fig. <ref>). Here, i, i+1, and i+2 take values in {1,2,3} and the sums should be interpreted mod 3: for instance, if i=2 then i+1=3 and i+2=1.

From (<ref>) we see that the change in the wedge angle enters into the energy not just quadratically but also linearly, therefore the variations of the wedge angle should be computed up to the second order approximation.
To keep the notation compact, we set

    _i:=(^ℓ+ℓ_i),    _0:=(^ℓ).

Let 

    _i+1   :=(^ℓ + ℓ_i+1 + (^ℓ+ℓ_i+1)) - (^ℓ +(^ℓ))
       =ℓ_i+1+(_i+1 -_0),

and

    _i+2   :=(^ℓ + ℓ_i+2 + (^ℓ+ℓ_i+2)) - (^ℓ +(^ℓ))
       =ℓ_i+2+(_i+2 -_0),

be the images of the edges parallel to _i+1 and _i+2 and starting at ^ℓ. Then, the angle ϑ_i=ϑ_i(^ℓ) is given by

    cos(ϑ_i)=_i+1·_i+2/|_i+1||_i+2|.
 


Calculations given in Appendix <ref> yield that

    ϑ_i= 2/3π+δϑ_i^(1)+δϑ_i^(2) + o(| |^2),

where δϑ_i^(1) and δϑ_i^(2) are the first order and the second order variation, respectively, of the wedge angle with respect to the reference angle 2/3π. Therefore, keeping up to second order terms one has that

    δϑ_i=δϑ_i^(1)+δϑ_i^(2).

It turns out that the first order variation takes the form
    δϑ_i^(1)(^ℓ)=-1/ℓ(_i+1-_0)·_i+1^⊥  + 1/ℓ(_i+2 -_0)·_i+2^⊥,

with  _i+1^⊥ defined by

    _i+1^⊥:=_i+2 +1/2_i+1/|_i+2 +1/2_i+1|=2/√(3)(_i+2 +1/2_i+1)      i=1, 2, 3,

that is, the unit vector orthogonal to _i+1, (cf. equation (19) in <cit.>).

Figure <ref> illustrates the geometrical meaning of formula (<ref>). 
In particular, 

    ∑_i=1^3δϑ_i^(1)(^ℓ)=0,

as it could have been deduced from geometrical considerations.

The second order variation is given by, see Appendix <ref>,

    δϑ_i^(2)(^ℓ)=-1/√(3)[   -1/2(δϑ_i^(1))^2+2/ℓ^2(𝐮_i+1 - 𝐮_0)·(𝐮_i+2 - 𝐮_0)-
       ((𝐮_i+1 - 𝐮_0)/ℓ^2·_i+1+(𝐮_i+2 - 𝐮_0)/ℓ^2·_i+2)×
       ×((𝐮_i+1 - 𝐮_0)/ℓ^2·_i+2+(𝐮_i+2 - 𝐮_0)/ℓ^2·_i+1 -√(3)/2δϑ_i^(1))
    
        -_i+1·_i+2/ℓ^4(|𝐮_i+1 - 𝐮_0|^2-21/ℓ^2( _i+1·(𝐮_i+1 - 𝐮_0) )^2
    
       +|𝐮_i+2 - 𝐮_0|^2-21/ℓ^2( _i+2·(𝐮_i+2 - 𝐮_0) )^2    ].

By algebric manipulation one finds that

    ∑_i=1^3δϑ_i(^ℓ)=∑_i=1^3δϑ_i^(2)(^ℓ)=-3√(3)/ℓ^2(1/3∑_i=1^3w(^ℓ+ℓ_i)-w(^ℓ))^2,

where w denotes the out-of-plane component of the displacement, that is

    w:=·_3,

where _3 is the unit vector perpendicular to the undeformed sheet.
Note that, by (<ref>), the ∑_iδϑ_i(^ℓ) is non-positive and hence the contribution of the self-stress to the strain energy is non-negative for τ_0<0, i.e., for δϑ_0>0, see (<ref>).





 §.§ Change of the dihedral angles


For each fixed node ^ℓ∈ L_2(ℓ)∩Ω and  for each edge parallel to _i
and starting at ^ℓ we need to define four types of dihedral angles __i^+(^ℓ), __i^-(^ℓ), __i_i+1(^ℓ) and __i_i+2(^ℓ):


    cos__i^+=(_i×_i+1)·(_i×_i^+)/|_i×_i+1||_i×_i^+|,
       cos__i^-=(_i+2×_i)·(_i^-×_i)/|_i+2×_i||_i^-×_i)|,
       cos__i_i+1=(_i×_i+1)·(_i^-×_i)/|_i×_i+1||_i^-×_i|,
       cos__i_i+2=(_i+2×_i)·(_i×_i^+)/|_i+2×_i||_i×_i^+|,

where	

    _i^+=   ^ℓ+ℓ_i-ℓ_i+2+(^ℓ+ℓ_i-ℓ_i+2)-( ^ℓ+ℓ_i+(^ℓ+ℓ_i) )
    
    =   -ℓ_i+2+_i^+-_i,         _i^+:=(^ℓ+ℓ_i-ℓ_i+2),
    _i^-=   ^ℓ+ℓ_i-ℓ_i+1+(^ℓ+ℓ_i-ℓ_i+1)-( ^ℓ+ℓ_i+(^ℓ+ℓ_i) )
    
    =   -ℓ_i+1+_i^--_i,         _i^-:=(^ℓ+ℓ_i-ℓ_i+1)

are the images of vectors ℓ_i^+ and ℓ_i^- (see Fig. <ref>, for i=1), parallel to _i+2 and _i+1 and starting at the image of the point ^ℓ+ℓ_i.


Also here, i, i+1, and i+2 take values in {1,2,3} and the sums should be interpreted mod 3: for instance, if i=3 then i+1=1 and i+2=2.

The  C-dihedral angle __i^+(^ℓ) is the angle corresponding to the C-dihedron with middle edge ℓ_i and oriented as _i^⊥,  while __i^-(^ℓ) is the angle corresponding to the C-dihedron oriented opposite to _i^⊥ (see Fig. <ref> for i=1). 
The Z-dihedral angle __i_i+1(^ℓ) 
corresponds to the Z-dihedron with middle edge ℓ_i and the other two edges parallel to _i+1
(see Fig. <ref> for i=1).


Then, recalling that δΘ=Θ, calculations in Appendix <ref> yield that

    δ__i^+(^ℓ)=2√(3)/3ℓ[2w(^ℓ)-w(^ℓ+ℓ_i+1)+w(^ℓ+ℓ_i-ℓ_i+2)-2w(^ℓ+ℓ_i)],
		
and

    δ__i_i+1(^ℓ)=
    2√(3)/3ℓ[w(^ℓ+ℓ_i-ℓ_i+1)-w(^ℓ+ℓ_i)+w(^ℓ+ℓ_i+1)-w(^ℓ)].
		
Analogous formulas hold for δ__i^- and  δ__i_i+2:

    δ__i^-(^ℓ)=-2√(3)/3ℓ[2w(^ℓ)-w(^ℓ+ℓ_i+2)+w(^ℓ+ℓ_i-ℓ_i+1)-2w(^ℓ+ℓ_i)],
       δ__i_i+2(^ℓ)=2√(3)/3ℓ[w(^ℓ+ℓ_i-ℓ_i+2)-w(^ℓ+ℓ_i)+w(^ℓ+ℓ_i+2)-w(^ℓ)].




§ SPLITTING OF THE ENERGY


The above calculations show that δϑ_i^(1) as well as δ l_i depend upon the in-plane components of , cf.  (<ref>) and (<ref>), while δϑ_i^(2), δ, and δ depend upon the out-of-plane component of , cf.  (<ref>), (<ref>), and (<ref>). This yields a splitting of the energy into membrane and bending parts 

    𝒰_ℓ=𝒰_ℓ^(m)+𝒰_ℓ^(b),   𝒰_ℓ^(b):=_ℓ^(s) +_ℓ^(d)
   
defined by

    𝒰_ℓ^(m):=1/2 ∑_ℰ k^l   (δ l)^2+1/2 ∑_𝒲 k^ϑ (δϑ^(1))^2 
       _ℓ^(s) :=τ_0∑_𝒲δϑ^(2), 
       _ℓ^(d) :=1/2 ∑_𝒵 k^  (δ)^2+1/2 ∑_𝒞 k^  (δ)^2,

where _ℓ^(s) is the self-energy (corresponding to the so-called cohesive energy in the literature) and _ℓ^(d) is the dihedral energy.
The analysis in a paper by Davini <cit.> applies here to the in-plane deformations, providing a continuum model of the graphene sheet within the framework of Γ-convergence theory. Hereafter, we concentrate on the out-of-plane deformations. 
With the notation introduced in Section <ref> we now write the bending energy more explicitly.
The self-energy can be written as

    𝒰^(s)_ℓ
    = ∑_^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ωτ_0 ∑_i = 
    	1^3δϑ_i^(2) (^ℓ),

where ∑_i =1^3δϑ_i^(2) (^ℓ) is given in (<ref>) in terms of the out-of-plane component of the displacement w.
We further split the dihedral energy  _ℓ^(d) in

    _ℓ^(d):=𝒰^𝒵_ℓ+𝒰^𝒞_ℓ,

where

    𝒰^𝒵_ℓ=1/2  k^𝒵∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (δ__i_i+2(^ℓ))^2+ (δ__i_i+1(^ℓ))^2

is the contribution of the Z-dihedra, and

    𝒰^𝒞_ℓ=1/2 k^𝒞 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (δ__i^+(^ℓ))^2+ (δ__i^-(^ℓ))^2

is the contribution of the C-dihedra. The Z- and C-dihedral angles appearing in (<ref>) and (<ref>) are given in terms of w in (<ref>)-(<ref>).

In the next section we deduce, by means of a formal analysis, a continuous version of the discrete bending energy 

    𝒰_ℓ^(b)=_ℓ^(s) +𝒰^𝒵_ℓ+𝒰^𝒞_ℓ,

from which we shall deduce
expressions for the sheet's bending stiffnesses. A rigorous analysis based on Γ-convergence theory will be done in a forthcoming paper <cit.>. 




§ THE CONTINUUM LIMIT


In this section we  find a continuous energy, defined over the domain Ω, that 
approximates the discrete bending energy 𝒰_ℓ^(b)  defined over the lattice (L_1(ℓ)∪ L_2(ℓ))∩Ω. This is achieved by letting the lattice size ℓ go to zero so that
(L_1(ℓ)∪ L_2(ℓ))∩Ω invades Ω.
With this in mind, in place of a function w:(L_1(ℓ)∪ L_2(ℓ))∩Ω→ℝ, we consider a twice continuously differentiable function w:Ω→ℝ.

Given two vectors  and , with ∂^2_w we denote the second partial derivative
of w in the directions /|| and /||, that is

    ∂^2_w=∇^2 w /||·/||,

where  ∇^2 w denotes the Hessian of w. Clearly, we also have

    ∂^2_w(x_0)=lim_ℓ→ 0w(x_0+ℓ+ℓ)-w(x_0+ℓ)-w(x_0+ℓ)+w(x_0)/ℓ^2 || ||.

The change of the Z-dihedra, see (<ref>)_2, can be rewritten, after setting

    _i:=_i-_i+2,

as

    δ__i_i+2(^ℓ)   =2√(3)/3ℓ[w(^ℓ+ℓ_i-ℓ_i+2)-w(^ℓ)-w(^ℓ+ℓ_i)+w(^ℓ+ℓ_i+2)]
       =2√(3)/3ℓ[w(^ℓ+ℓ_i)-w(^ℓ)-w(^ℓ+ℓ_i+ℓ_i+2)+w(^ℓ+ℓ_i+2)]
       =2√(3)/3ℓ [- ∂^2__i_i+2w(^ℓ) ℓ^2 |_i| |_i+2|+o(ℓ^2)]
       =-2 ℓ∂^2__i_i+2w(^ℓ) +o(ℓ),

where the third equality follows from (<ref>). Similarly, setting

    _i:=_i-_i+1,

we have that

    δ__i_i+1(^ℓ)   =
    2√(3)/3ℓ[w(^ℓ+ℓ_i-ℓ_i+1)-w(^ℓ+ℓ_i)+w(^ℓ+ℓ_i+1)-w(^ℓ)]
       =2√(3)/3ℓ[w(^ℓ+ℓ_i)-w(^ℓ+ℓ_i+ℓ_i+1)+w(^ℓ+ℓ_i+1)-w(^ℓ)]
       =-2 ℓ∂^2__i_i+1w(^ℓ) +o(ℓ).

Taking (<ref>) into account, we may rewrite the vectors _i and _i, defined in (<ref>) and (<ref>), in terms of _i, for instance _1=_1 and _1=-_3,
and then rewrite the Z-dihedral energy, see (<ref>) and rewritten below for the reader convenience, as

    𝒰^𝒵_ℓ   =1/2  k^𝒵∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (δ__i_i+2(^ℓ))^2+ (δ__i_i+1(^ℓ))^2
    
       =1/24ℓ^2 k^𝒵∑_^ℓ∈ L_2(ℓ)∩Ω (∂^2__1_3w(^ℓ))^2+ (∂^2__1_1w(^ℓ))^2+(∂^2__2_2w(^ℓ))^2
       +(∂^2__2_3w(^ℓ))^2+(∂^2__3_1w(^ℓ))^2+(∂^2__3_2w(^ℓ))^2+o(ℓ^2),
       =1/28√(3)/9 k^𝒵∑_^ℓ∈ L_2(ℓ)∩Ω(∂^2__1_3w(^ℓ))^2+ (∂^2__1_1w(^ℓ))^2+(∂^2__2_2w(^ℓ))^2
       +(∂^2__2_3w(^ℓ))^2+(∂^2__3_1w(^ℓ))^2+(∂^2__3_2w(^ℓ))^2)|E^ℓ(^ℓ)|+o(ℓ^2),

where |E^ℓ(^ℓ)|=ℓ^23√(3)/2 is the area of the hexagon 
E^ℓ(^ℓ) of side ℓ centred at ^ℓ (see Fig. <ref>).





Let χ_E^ℓ(^ℓ)() be the characteristic function of E^ℓ(^ℓ), i.e., it is equal to 1 if ∈  
E^ℓ(^ℓ) and 0 otherwise, and let 

    W^_ℓ():=∑_^ℓ∈ L_2(ℓ)∩Ω   ((∂^2__1_3w(^ℓ))^2+ (∂^2__1_1w(^ℓ))^2+(∂^2__2_2w(^ℓ))^2
       +(∂^2__2_3w(^ℓ))^2+(∂^2__3_1w(^ℓ))^2+(∂^2__3_2w(^ℓ))^2)χ_E^ℓ(^ℓ)().

Then, we may simply write

    𝒰^𝒵_ℓ=1/28√(3)/9 k^𝒵∫_Ω W^_ℓ() d+o(ℓ^2),

and since W^_ℓ converges, as ℓ goes to zero, to

    (∂^2__1_3w)^2+ (∂^2__1_1w)^2+(∂^2__2_2w)^2+(∂^2__2_3w)^2+(∂^2__3_1w)^2+(∂^2__3_2w)^2,

we deduce that

    lim_ℓ→ 0𝒰^𝒵_ℓ=1/28√(3)/9 k^𝒵∫_Ω   (∂^2__1_3w)^2+ (∂^2__1_1w)^2
    +(∂^2__2_2w)^2
       +(∂^2__2_3w)^2
    +(∂^2__3_1w)^2+(∂^2__3_2w)^2 d=:𝒰^𝒵_0(w).

The functional 𝒰^𝒵_0, defined in (<ref>), is the continuum limit of the Z-dihedral energy.



Working in a similar manner, we find the continuum limit of the C-dihedral energy:

    lim_ℓ→ 0𝒰^𝒞_ℓ=1/216√(3)/9 k^𝒞∫_Ω∑_i=1^3 (∂^2__i_i^⊥w)^2 d
    =:𝒰^𝒞_0(w),

and the continuum limit of the self-energy:

    lim_ℓ→ 0𝒰^(s)_ℓ=-4/9τ_0 ∫_Ω ( ∂^2__1 _1w+∂^2__2 _2w+∂^2__1 _2w)^2  d
    =:𝒰^(s)_0(w).

Detailed calculations leading to (<ref>) and (<ref>) are found in Appendix <ref>.


The total bending limit energy is therefore

    𝒰^(b)_0(w):=𝒰^𝒵_0(w)+𝒰^𝒞_0(w)+𝒰^(s)_0(w).





§ THE EQUIVALENT PLATE EQUATION

In this section we rewrite the limit energies in a more amenable form.

We start by manipulating the limit C-dihedral energy. We first note that

    ∂^2__2_2^⊥w   =∇^2 w _2·_2^⊥=∇^2 w (-1/2_1+√(3)/2_1^⊥)· (-√(3)/2_1-1/2_1^⊥)
       =√(3)/4∂^2__1_1w-√(3)/4∂^2__1^⊥_1^⊥w
    -1/2∂^2__1_1^⊥w,

and similarly

    ∂^2__3_3^⊥w=
    -√(3)/4∂^2__1_1w+√(3)/4∂^2__1^⊥_1^⊥w
    -1/2∂^2__1_1^⊥w,

from which we find that 

    ∑_i=1^3 (∂^2__i_i^⊥w)^2   =(∂^2__1_1^⊥w)^2+
    (√(3)/4∂^2__1_1w-√(3)/4∂^2__1^⊥_1^⊥w
    -1/2∂^2__1_1^⊥w)^2
         +
    (-√(3)/4∂^2__1_1w+√(3)/4∂^2__1^⊥_1^⊥w
    -1/2∂^2__1_1^⊥w)^2
       =3/2 (∂^2__1_1^⊥w)^2+ 3/8(∂^2__1_1w-∂^2__1^⊥_1^⊥w)^2
       =3/2 (∂^2__1_1^⊥w)^2+ 3/8(∂^2__1_1w+∂^2__1^⊥_1^⊥w)^2-3/2∂^2__1_1w∂^2__1^⊥_1^⊥w
       =3/8 (Δ w)^2- 3/2∇^2 w,

where Δ w denotes the Laplacian of w.
Hence, the C-dihedral energy defined in (<ref>) rewrites as

    𝒰^𝒞_0(w)   =1/216√(3)/9 k^𝒞∫_Ω∑_i=1^3 (∂^2__i_i^⊥w)^2 d=
    1/22√(3)/3 k^𝒞∫_Ω (Δ w)^2- 4 ∇^2 w d.





We now tackle the Z-dihedral energy. Recalling (<ref>), by a calculation similar to that carried on in (<ref>) we find

    ∂^2__1_3w=
    -√(3)/2∂^2__1_1w-1/2∂^2__1_1^⊥w     ∂^2__1_1w=
    +√(3)/2∂^2__1_1w-1/2∂^2__1_1^⊥w,

where _i^⊥=_3×_i, and from these equations we deduce that

    (∂^2__1_1w)^2+(∂^2__1_3w)^2=3/2 (∂^2__1_1w)^2+
    1/2 (∂^2__1_1^⊥w)^2.

Similar identities hold for _2 and _3. Thence, we find that  the Z-dihedral energy
takes the form

    𝒰^𝒵_0(w)   =1/28√(3)/9 k^𝒵∫_Ω (∂^2__1_3w)^2+ (∂^2__1_1w)^2
    +(∂^2__2_2w)^2
       +(∂^2__2_3w)^2
    +(∂^2__3_1w)^2+(∂^2__3_2w)^2 d
       =1/28√(3)/9 k^𝒵∫_Ω3/2∑_i=1^3(∂^2__i_iw)^2+
    1/2∑_i=1^3(∂^2__i_i^⊥w)^2  d.

The second sum is equal, as it can be checked, to the last line of (<ref>), and with a similar calculation we also find that

    ∑_i=1^3(∂^2__i_iw)^2=9/8 (Δ w)^2- 3/2∇^2 w,

and hence

    𝒰^𝒵_0(w)
       =1/28√(3)/9 k^𝒵∫_Ω3/2(9/8 (Δ w)^2- 3/2∇^2 w)+
    1/2(3/8 (Δ w)^2- 3/2∇^2 w)  d
       =1/28√(3)/9 k^𝒵∫_Ω15/8 (Δ w)^2- 3 ∇^2 w d
       =1/25√(3)/3 k^𝒵∫_Ω  (Δ w)^2- 8/5∇^2 w d.



We now deal with the self-energy. Again with a calculation similar to that carried on in (<ref>) we find

    ∂^2__1_2w
       =-1/2∂^2__1_1w+√(3)/2∂^2__1_1^⊥w
    
    ∂^2__2_2w   =1/4∂^2__1_1w+3/4∂^2__1^⊥_1^⊥w
    -√(3)/2∂^2__1_1^⊥w,

and hence

    𝒰^(s)_0(w)   =-4/9τ_0 ∫_Ω ( ∂^2__1 _1w+∂^2__2 _2w+∂^2__1 _2w)^2  d
       =-4/9τ_0 ∫_Ω ( 3/4∂^2__1 _1w+ 3/4∂^2__1^⊥_1^⊥w)^2  d
       =-1/2τ_0/2∫_Ω ( Δ w)^2  d.


By summing (<ref>), (<ref>), and (<ref>), we find the total bending energy, defined in (<ref>):

    𝒰^(b)_0(w)
    =   1/2∫_Ω(5√(3)/3 k^𝒵+2√(3)/3 k^𝒞-τ_0/2)( Δ w)^2
       +(- 8/55√(3)/3 k^𝒵-42√(3)/3 k^𝒞)∇^2 w d,
       =1/2∫_Ω ( Δ w)^2+_G∇^2 w d,

where 

    := 5√(3)/3 k^𝒵+2√(3)/3 k^𝒞-τ_0/2,     _G:=- 8/55√(3)/3 k^𝒵-42√(3)/3 k^𝒞.

are the bending and the Gaussian stiffnesses, respectively. _G is called Gaussian because it multiplies the Gaussian curvature ∇^2w, while Δ w is twice the mean curvature. The analytical expression for the bending stiffness (<ref>)_1 coincides with the one deduced in <cit.>, within a discrete mechanical framework, if one assumes that k^𝒵≡ k^𝒞. It clearly shows that the origin of the bending stiffness is twofold: a part depends on the dihedral contribution, and a part on the self-stress. Quite surprisingly, the self-stress has no role in the Gaussian stiffness. It is worth noticing that in the above approach there is no need of introducing any questionable effective thickness parameter.





§ NUMERICAL RESULTS

In this section, we adopt the 2nd-generation Brenner potential <cit.> to obtain quantitative results for the continuum material parameters deduced in Sec. <ref>.

The 2nd-generation REBO potentials developed for hydrocarbons by Brenner et al. in <cit.> accommodate  up to third-nearest-neighbor interactions through a bond-order function depending, in particular, on dihedral angles.  Following Appendix B of <cit.>, we here give a short account of the form of this potential.

The  binding energy V of an atomic aggregate is given as a sum over nearest neighbors:


    V=∑_i∑_j<i V_ij ;


the interatomic potential V_ij is given by 


    V_ij=V_R(r_ij)+b_ijV_A(r_ij),


where the individual effects of the repulsion and attraction functions V_R(r_ij) and V_A(r_ij), which model pair-wise interactions of  atoms i and j depending on their distance r_ij, are modulated by the bond-order function b_ij. The repulsion and attraction functions have the following forms:


    V_A(r)   =-f^C(r)∑_n=1^3B_n e^-β_n r ,
    
    V_R(r)   =f^C(r)( 1 + Q/r) A e^-α r ,

where f^C(r) is a cutoff function limiting the range of covalent interactions, and where Q, A, B_n, α, and β, are parameters to be chosen fit to some material-specific dataset. The remaining ingredient in (<ref>) is the bond-order function:


    b_ij=1/2(b_ij^σ-π+b_ji^σ-π)+b_ij^π ,


where apexes σ and π refer to two types of bonds: the strong covalent σ-bonds between atoms in one and the same given plane, and the π-bonds responsible for interlayer interactions, which are perpendicular to the plane of σ-bonds. 
The role of function b_ij^σ-π is to account for the local coordination of, and the bond angles relative to, atoms i and j; its form is:


    b_ij^σ-π=(1+∑_k≠ i,j f_ik^C(r_ik)G(cosθ_ijk)  e^λ_ijk+P_ij(N_i^C,N_i^H)  )^-1/2 .


Here, for each fixed pair of indices (i,j), (a) the cutoff function f_ik^C limits the interactions of atom i to those with its nearest neighbors; (b) λ_ijk is a string of parameters designed to prevent attraction in some specific situations; (c) function P_ij depends on N_i^C and N_i^H, the numbers of C and H atoms that are nearest neighbors of atom i; it is meant to adjust the bond-order function 	according to the environment of the C atoms in one or another molecule; (d) for solid-state carbon, the values of  both the string λ_ijk and the function P_ij are taken null; (e)  function G modulates the contribution of each nearest neighbour of atom i in terms of the cosine of the angle between the ij and ik bonds; its analytic form is given by three sixth-order polynomial splines.	
Function b_ij^π is given a split representation:


    b_ij^π=Π_ij^RC+b_ij^DH,

where the first addendum Π_ij^RC depends on whether the bond between atoms i and j has a radical character and on whether it is part of a conjugated system, while the second addendum b_ij^DH depends on dihedral angles and has the following form:

    b_ij^DH=T_ij(N_i^t,N_j^t,N_ij^ conj)(∑_k(≠ i,j)∑_k(≠ i,j)( 1-cos^2Θ_ijkl)f_ik^C(r_ik)f_jl^C(r_jl)  ) ,

where function T_ij is a tricubic spline depending on N_i^t=N_i^C+N_i^H, N_j^t, and N_ij^ conj, a function  of local conjugation, and the dihedral angle is defined as

    cosΘ_ijkl=_jik·_ijl,   _jik=_ji×_ik/|_ji×_ik|,  _jil=_ij×_il/|_ij×_jl|.


The values of the constant k^ and k^ can be deduced by deriving twice of the potential, and computing the result in the ground state (GS): r_ij=ℓ, θ_ijk=2/3π, Θ_ijkl=0. In particular, we find:

    k^Θ:=k^=k^=∂^2_Θ_ijklV_ij|_GS=2TV_A(ℓ),

where T is the value of T_ij in the GS. 


	With this notation the bending stiffness becomes:
	
    =7√(3)/3k^Θ-τ_0/2.

	This expression coincides with that given in [28]:
	
    𝒟=V_A(r_0)/2(  (b_0^σ-π)'-14 T_0/√(3)),

	after noticing that 
	
    V_A(r_0)(b_0^σ-π)'≡-τ_0,    -V_A(r_0)(  7 T_0/√(3))≡ 2TV_A(ℓ) 7/√(3)= 7/√(3)k^Θ.

	
	Instead, in references [4] and [21] the dihedral energies are not contemplated and the bending stiffness found, up to notational differences, coincides with   ours after setting k^Θ=0.




With the values reported in <cit.>, we get:

    k^Θ=0.0282 nN nm=0.1764  eV.

From <cit.>, we take the value of the selfstress τ_0:

    τ_0=-0.2209 nN nm=-1.3787 eV.

With (<ref>) and (<ref>), we obtain:

    =7√(3)/3k^Θ-τ_0/2=0.2247  nN nm=1.4022  eV,
       _G=-16√(3)/3k^Θ=-0.2610 nN nm=-1.6293 eV.

The value  of  is in complete agreement with the literature <cit.>; from (<ref>) and (<ref>) it is possible to check that  the contribution of the self-stress and the dihedral stiffness amounts to about  49.16% and 50.84% of the total, respectively. Neither analytical evaluation of _G, nor MD computations, have been proposed so far. The value we obtain is in good agreement with the value of -1.52 eV, reported in <cit.> and  determined by means of DFT.




§ CONCLUSIONS


Starting form a discrete model inferred from MD, we have deduced a continuum theory describing the bending behavior of a graphene sheet. Atomic interactions have been modeled by exploiting the main features of the 2nd-generation Brenner potential and adopting a quadratic approximation of the energy. The deduced continuum limit   fully describes the bending behavior of graphene. To our knowledge, it is the first time that an analytical expression of the Gaussian stiffness is given  and an explanation of its origins at the atomistic scale is provided. We also derived a quantitative evaluation of the related constitutive parameters.



	

§ ACKNOWLEDGMENTS
	
	AF acknowledges the financial support of Sapienza University of Rome (Progetto d'Ateneo 2016 — “Multiscale Mechanics of 2D Materials: Modeling and Applications”). 
	





§ APPENDIX

In order to compute the strain measures for small changes of configuration of the graphene foil, we write the displacements of the nodes in the form

    u (^ℓ)= ξ(^ℓ),

where ξ is a positive scalar measuring smallness and := u/ξ stands for the displacement distribution normalized accordingly. 






 §.§ Change of the bond angle


Let us define the bond angle as

    cos(ϑ_i(ξ))=(𝐦(ξ)·𝐧(ξ)/|𝐦(ξ)||𝐧(ξ)|) ,
 
with

    𝐦(ξ):=ℓ_i+1+ξ (_i+1 - _0)    𝐧(ξ):=ℓ_i+2+ξ (_i+2 - _0),

where we have set

    _0:=(^ℓ),   _i:=(^ℓ+ℓ_i).

Then, from Taylor's expansion we get

    ϑ_i(ξ)= ϑ_i(0)+ϑ_i^'(0) ξ+1/2 ϑ_i^''(0) ξ^2 +O(ξ^3),

where the various terms can be calculated by successive differentiations of Eq. (<ref>). Thus,

    -sin(ϑ_i(0))    ϑ_i^'(0)=.(𝐦(ξ)·𝐧(ξ)/|𝐦(ξ)||𝐧(ξ)|)^'|_ξ=0
       =𝐦'(ξ)·𝐧(ξ)+𝐦(ξ)·𝐧'(ξ)/|𝐦(ξ)||𝐧(ξ)|-𝐦(ξ)·𝐧(ξ)/|𝐦(ξ)||𝐧(ξ)|.(𝐦(ξ)·𝐦'(ξ)/|𝐦(ξ)|^2 +𝐧(ξ)·𝐧'(ξ)/|𝐧(ξ)|^2)|_ξ=0,

which yields

    ϑ_i^'(0)=-_i+2 +1/2_i+1/|_i+2 +1/2_i+1|·(_i+1 - _0) - _i+1 +1/2_i+2/|_i+1 +1/2_i+2|·(_i+2 - _0),

where we take into account that sin(ϑ_i(0))=√(3)/2, |𝐦(0)|=|𝐧(0)|=ℓ, |_i+1 +1/2_i+2|=|_i+2 +1/2_i+1| = √(3)/2 and _i·_i+1=-1/2.

Moreover, by differentiating Eq. (<ref>) twice, we get

    -cos(ϑ_i(0)) ϑ^'_i(0)^2-sin(ϑ_i(0)) ϑ_i^''(0)=. (𝐦(ξ)·𝐧(ξ)/|𝐦(ξ)||𝐧(ξ)|)^'' |_ξ=0,

which gives

    ϑ_i^''(0)= -1/sinϑ_i(0)(cosϑ_i(0)ϑ^'_i(0)^2+.(𝐦(ξ)·𝐧(ξ)/|𝐦(ξ)||𝐧(ξ)|)”|_ξ=0).
 
Computations yield that

    (𝐦·𝐧/|𝐦||𝐧|)”   =𝐦”/|𝐦|·𝐧/|𝐧|+2 𝐦'/|𝐦|·𝐧'/|𝐧|+𝐦/|𝐦|·𝐧”/|𝐧|
       -[𝐦'/|𝐦|·𝐧/|𝐧|+𝐦/|𝐦|·𝐧'/|𝐧|+(𝐦/|𝐦|·𝐧/|𝐧|)'  ] (𝐦'/|𝐦|·𝐦/|𝐦|+𝐧'/|𝐧|·𝐧/|𝐧|)
       -𝐦/|𝐦|·𝐧/|𝐧|[ (|𝐦'|/|𝐦|)^2+𝐦/|𝐦|·𝐦”/|𝐦|-2(𝐦'/|𝐦|·𝐦/|𝐦|)^2.
       . +
    (|𝐧'|/|𝐧|)^2+𝐧/|𝐧|·𝐧”/|𝐧|-2(𝐧'/|𝐧|·𝐧/|𝐧|)^2
    ] .

Since 𝐦”(0)=𝐧”(0)=0 and (𝐦·𝐧/|𝐦||𝐧|)'|_ξ=0=-sinϑ_i(0)ϑ_i^'(0), we finally have that

    .(𝐦·𝐧/|𝐦||𝐧|)”|_ξ=0   =2/ℓ^2(_i+1 - _0)·(_i+2 - _0)-((_i+1 - _0)/ℓ·_i+1+(_i+2 - _0)/ℓ·_i+2)×
       ×((_i+1 - _0)/ℓ·_i+2+(_i+2 - _0)/ℓ·_i+1 -sinϑ_i(0)ϑ_i^'(0)  )
        -(_i+1·_i+2)(1/ℓ^2|_i+1 - _0|^2-21/ℓ^2( _i+1·(_i+1 - _0) )^2
    +1/ℓ^2|_i+2 - _0|^2 .
       .-21/ℓ^2( _i+2·(_i+2 - _0) )^2 
    ).





All in all,  we have that:

    ϑ_i”(0)=-2/√(3)[-1/2ϑ_i'(0) ^2+2/ℓ^2(_i+1 - _0)·(_i+2 - _0)-
       ((_i+1 - _0)/ℓ·_i+1+(_i+2 - _0)/ℓ·_i+2)×
       ×((_i+1 - _0)/ℓ·_i+2+(_i+2 - 𝔲_0)/ℓ·_i+1 -√(3)/2ϑ_i^'(0)  )
        -_i+1·_i+2/ℓ^2(|_i+1 - _0|^2-2( _i+1·(_i+1 - _0) )^2
    +|_i+2 - _0|^2 
       -2( _i+2·(_i+2 - _0) )^2    ] .

Recalling that δϑ_i=ϑ_i^'(0) ξ+1/2 ϑ_i^''(0) ξ^2 +O(ξ^3) and that 𝐮=ξ, we get:

    δϑ_i=δϑ_i^(1)+δϑ_i^(2)+O(ξ^3),

with

    δϑ_i^(1)=ϑ_i^'(0)ξ=-1/ℓ_i+2 +1/2_i+1/|_i+2 +1/2_i+1|·(𝐮_i+1 - 𝐮_0) - 1/ℓ_i+1 +1/2_i+2/|_i+1 +1/2_i+2|·(𝐮_i+2 - 𝐮_0).

and

    δϑ_i^(2)=1/2ϑ_i^''(0)ξ^2= -1/√(3)[-1/2(ξϑ_i'(0)) ^2+2/ℓ^2(𝐮_i+1 - 𝐮_0)·(𝐮_i+2 - 𝐮_0)-
       ((𝐮_i+1 - 𝐮_0)/ℓ^2·_i+1+(𝐮_i+2 - 𝐮_0)/ℓ^2·_i+2)×
       ×((𝐮_i+1 - 𝐮_0)/ℓ^2·_i+2+(𝐮_i+2 - 𝐮_0)/ℓ^2·_i+1 -√(3)/2ξϑ_i^'(0)  )
    
        -_i+1·_i+2/ℓ^4(|𝐮_i+1 - 𝐮_0|^2-21/ℓ^2( _i+1·(𝐮_i+1 - 𝐮_0) )^2
    +|𝐮_i+2 - 𝐮_0|^2 
       -21/ℓ^2( _i+2·(𝐮_i+2 - 𝐮_0) )^2    ].

We write (<ref>) in the simpler form

    δϑ_i^(1)=-1/ℓ (𝐮_i+1 - 𝐮_0)·^⊥_i+1 + 1/ℓ (𝐮_i+2 - 𝐮_0)·^⊥_i+2,

where the unit vectors ^⊥_i+1's are defined by (<ref>). 



 §.§ Change of the dihedral  angle

To fix the ideas, we focus on the  dihedral angles δ__1^+, and  δ__1_2, sketched in Fig. <ref>; the other strains can be obtained in analogous manner.  


The first order approximation of the dihedral angle is all we need to evaluate the corresponding energy contribution. 

Let us denote by 

    _1=ℓ_1+(_1-_0)ξ,   _2=ℓ_2+(_2-_0)ξ,   _4=ℓ_4+(_4-_1)ξ

the edge vectors after the deformation. We have that

    |sin( __1^+(ξ))|=|(_1×_2)×(_1×_4)|/|_1×_2||_1×_4|;

it is easy to see that

    _1×_2=√(3)/2ℓ^2_3+ ( ℓ_1×(_2-_0) -ℓ_2×(_1-_0))ξ+O(ξ^2),
       _1×_4=√(3)/2ℓ^2_3+( ℓ_1×(_4-_1) -ℓ_4×(_1-_0))ξ+O(ξ^2),

whence

    (_1×_2)×(_1×_4)=√(3)/2_3   ×( ℓ_1×(_4-_1)-ℓ_4×(_1-_0) +
       +ℓ_2×(_1-_0)-ℓ_1×(_2-_0))ξ+O(ξ^2).

On recalling the  identity

    ×(×)=(·)-(·),

we obtain:

    (_1×_2)×(_1×_4)=√(3)/2ℓ^2(2𝔴_1-𝔴_4+𝔴_2-2𝔴_0 )ℓ_1ξ+O(ξ^2),

where we have used the fact that _4-_2(=-_3-_2)=_1 and set 𝔴=·_3.
Now,

    |_1×_2|=|_1×_2|=√(3)/2ℓ^2+O(ξ^2),

whence

    |sin( __1^+(ξ))|=|δ__1^+|ξ+O(ξ^2)=2√(3)/3ℓ|2𝔴_1-𝔴_4+𝔴_2-2𝔴_0  |ξ+O(ξ^2).

Thus, on recalling that w=𝔴ξ, we conclude that

    δ__1^+=2√(3)/3ℓ(2w_1-w_4+w_2-2w_0  ).

For a generic C-dihedral angle centered in _i, we get

    δ__i^+(^ℓ)=2√(3)/3ℓ[2w(^ℓ)-w(^ℓ+ℓ_i+1)+w(^ℓ+ℓ_i-ℓ_i+2)-2w(^ℓ+ℓ_i)],
       δ__i^-(^ℓ)=-2√(3)/3ℓ[2w(^ℓ)-w(^ℓ+ℓ_i+2)+w(^ℓ+ℓ_i-ℓ_i+1)-2w(^ℓ+ℓ_i)].
		
For the Z-dihedral angle __1_2, we introduce the vector 

    _5=ℓ_5+(_5-_1)ξ,

image of _5 under the deformation. We have that 

    |sin__1_2|=|(_1×_2)×(_5×_1)|/|_1×_2||_5×_1|,

and

    _5×_1=√(3)/2ℓ^2_3+( ℓ_5×(_1-_0)-ℓ_1×(_5-_1) )ξ+O(ξ^2),

whence

    (_1×_2)×(_5×_1)=√(3)/2ℓ^2_3   ×( ℓ_5×(_1-_0)-ℓ_1×(_5-_1)
       +ℓ_2×(_1-_0)-ℓ_1×(_2-_0) )ξ+O(ξ^2).

Again, on making use of (<ref>) and recalling that _5=-_2, we obtain:

    (_1×_2)×(_5×_1)=√(3)/2ℓ^2(𝔴_5-𝔴_1+𝔴_2-𝔴_0)ξ+O(ξ^2).

On noticing that |_5×_1|=√(3)/2ℓ^2+O(ξ^2), we get

    |sin__1_2|=|δ__1_2|ξ+O(ξ^2)=2√(3)/3ℓ|𝔴_5-𝔴_1+𝔴_2-𝔴_0  |ξ+O(ξ^2),

whence

    δ__1_2=2√(3)/3ℓ(w_5-w_1+w_2-w_0  ).

For a generic Z-dihedral angle centered in _i, we get

    δ__i_i+1(^ℓ)=2√(3)/3ℓ[w(^ℓ+ℓ_i-ℓ_i+1)-w(^ℓ+ℓ_i)+w(^ℓ+ℓ_i+1)-w(^ℓ)],
       δ__i_i+2(^ℓ)=2√(3)/3ℓ[w(^ℓ+ℓ_i-ℓ_i+2)-w(^ℓ+ℓ_i)+w(^ℓ+ℓ_i+2)-w(^ℓ)].





 §.§ Deduction of the continuum limits 𝒰^𝒞_0 and 𝒰^(s)_0


In this Appendix we give a justification of (<ref>) and (<ref>) following the lines outlined in Section <ref> for the derivation of the continuum limit of the Z-dihedral energy.


As in Section <ref>,  in place of a function w:(L_1(ℓ)∪ L_2(ℓ))∩Ω→ℝ, we consider a twice continuously differentiable function w:Ω→ℝ.
We first consider the contribution due to the C-dihedra.
Momentarily, to keep the notation compact, we set

    :=∇ w(^ℓ),     :=∇^2 w(^ℓ).

By Taylor's expansion we have

    -w(^ℓ+ℓ_i+1)   =-w(^ℓ)-ℓ·_i+1-1/2ℓ^2 _i+1·_i+1+o(ℓ^2)
    
    w(^ℓ+ℓ_i-ℓ_i+2)   =w(^ℓ)+ℓ· (_i-_i+2)+1/2ℓ^2  (_i-_i+2)· (_i-_i+2)+o(ℓ^2)
    
    -2w(^ℓ+ℓ_i)   =-2w(^ℓ)-2ℓ·_i- ℓ^2 _i·_i+o(ℓ^2).

Thence, since _1+_2+_3=0, we find from (<ref>) that

    δ__i^+(^ℓ)   =2√(3)/3ℓ[2w(^ℓ)-w(^ℓ+ℓ_i+1)+w(^ℓ+ℓ_i-ℓ_i+2)-2w(^ℓ+ℓ_i)]
       =√(3)ℓ/3[-_i+1·_i+1+ (_i-_i+2)· (_i-_i+2)-2_i·_i)]+o(ℓ)

and, by substituting _i+2=-(_i+_i+1), we eventually get

    δ__i^+(^ℓ)   =√(3)ℓ/3[-_i+1·_i+1+ (2_i+_i+1)· (2_i+_i+1)-2_i·_i]+o(ℓ)
       =4√(3)ℓ/3[_i· (1/2_i+ _i+1)]+o(ℓ)
       =2ℓ_i·_i^⊥+o(ℓ)=2ℓ∂^2__i_i^⊥w(^ℓ)+o(ℓ).


where _i^⊥ is defined in (<ref>). Similarly, we also find that

    δ__i^-(^ℓ)=2ℓ∂^2__i_i^⊥w(^ℓ)+o(ℓ).

With the same steps taken in the study of the Z-dihedral energy we now derive the limit of the C-dihedral energy.
With the expressions of the change of the C-dihedra (<ref>) and (<ref>), we can rewrite the C-dihedral energy,
see (<ref>), as

    𝒰^𝒞_ℓ   =1/2 k^𝒞 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (δ__i^+(^ℓ))^2+ (δ__i^-(^ℓ))^2
    
       =1/2 8ℓ^2 k^𝒞 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (∂^2__i_i^⊥w(^ℓ))^2+o(ℓ^2)
    
       =1/216√(3)/9 k^𝒞 ∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (∂^2__i_i^⊥w(^ℓ))^2|E^ℓ(^ℓ)|+o(ℓ^2)
    
       =1/216√(3)/9 k^𝒞∫_Ω W^_ℓ() d+o(ℓ^2),

where the function

    W^_ℓ():=∑_^ℓ∈ L_2(ℓ)∩Ω∑_i=1^3 (∂^2__i_i^⊥w(^ℓ))^2 χ_E^ℓ(^ℓ)()

converges to ∑_i=1^3 (∂^2__i_i^⊥w)^2, as ℓ goes to zero.
Thus, 

    lim_ℓ→ 0𝒰^𝒞_ℓ=1/216√(3)/9 k^𝒞∫_Ω∑_i=1^3 (∂^2__i_i^⊥w)^2 d
    =:𝒰^𝒞_0(w),

which is (<ref>).

We now compute the limit of the self-energy.
By Taylor's expansion, with the notation introduced in (<ref>), and taking into account that
_1+_2+_3=0, we find

    ∑_iδϑ_i(^ℓ)   =-3√(3)/ℓ^2[1/3∑_i=1^3w(^ℓ+ℓ_i)-w(^ℓ)]^2
       =-3√(3)/ℓ^2ℓ^4/36 ( _1·_1+_2·_2+_3·_3)^2+o(ℓ^2)
       =-3√(3)/ℓ^2ℓ^4/9( _1·_1+_2·_2+_1·_2)^2+o(ℓ^2)
       =-√(3)/3ℓ^2( ∂^2__1 _1w(^ℓ)+∂^2__2 _2w(^ℓ)+∂^2__1 _2w(^ℓ))^2+o(ℓ^2).

With this expression the self-energy (<ref>) takes the form

    𝒰^(s)_ℓ   = ∑_^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ωτ_0 ∑_i = 
    	1^3δϑ_i (^ℓ) 
       = -√(3)/3τ_0 ℓ^2∑_^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω( ∂^2__1 _1w(^ℓ)+∂^2__2 _2w(^ℓ)+∂^2__1 _2w(^ℓ))^2+o(ℓ^2).

Since the sum is over the points of both lattices, whereas in the previous cases the sum was only over the nodes of L_2(ℓ), we cannot use the hexagons E^ℓ(^ℓ) earlier introduced. Let T^ℓ(^ℓ) be the triangle 
centered at ^ℓ of side √(3)ℓ as depicted in Figure <ref>.


Let

    W^ϑ_ℓ():=∑_^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω(∂^2__1 _1w(^ℓ)+∂^2__2 _2w(^ℓ)+∂^2__1 _2w(^ℓ)
    )^2 χ_T^ℓ(^ℓ)()

and note that the area of T^ℓ(^ℓ) is |T^ℓ(^ℓ)|=3√(3)/4ℓ^2.
The self-energy rewrites as

    𝒰^(s)_ℓ   = -4/9τ_0 ∑_^ℓ∈ (L_1(ℓ)∪ L_2(ℓ))∩Ω( ∂^2__1 _1w(^ℓ)+∂^2__2 _2w(^ℓ)+∂^2__1 _2w(^ℓ))^2|T^ℓ(^ℓ)| 
         +o(ℓ^2) 
       = -4/9τ_0 ∫_Ω W^ϑ_ℓ()  d+o(ℓ^2).

Since W^ϑ_ℓ converges to ( ∂^2__1 _1w+∂^2__2 _2w+∂^2__1 _2w)^2 as ℓ goes to zero, we find

    lim_ℓ→ 0𝒰^(s)_ℓ=-4/9τ_0 ∫_Ω ( ∂^2__1 _1w+∂^2__2 _2w+∂^2__1 _2w)^2  d
    =:𝒰^(s)_0(w),

which is (<ref>).

10
	
	Akiwande_2016
	D. Akinwande, C.J. Brennan, J.S. Bunch, P. Egberts, J.R. Felts, H. Gao,
	R. Huang, J.-S. Kim, T. Li, Y. Li, K.M. Liechti, N. Lu, H.S. Park, E.J. Reed,
	P. Wang, B.I. Yakobson, T. Zhang, Y.-W. Zhang, Y. Zhou, and Zhu Y.
	A review on mechanics and mechanical properties of 2d materials -
	graphene and beyond.
	Preprint at arXiv:1609.07187.
	
	Alessi_2016
	R. Alessi, A. Favata, and A. Micheletti.
	Pressurized CNTs under tension: A finite-deformation
	lattice model.
	Compos. Part B Eng., in press,
	http://dx.doi.org/10.1016/j.compositesb.2016.10.006.
	
	Arroyo_2002
	M. Arroyo and T. Belytschko.
	An atomistic-based finite deformation membrane for single layer
	crystalline films.
	J. Mech. Phys. Solids, 50(9):1941 – 1977, 2002.
	
	Arroyo2004
	M. Arroyo and T. Belytschko.
	Finite crystal elasticity of carbon nanotubes based on the
	exponential cauchy-born rule.
	Phys. Rev. B, 69:115415, 2004.
	
	Bajaj_2013
	C. Bajaj, A. Favata, and P. Podio-Guidugli.
	On a nanoscopically-informed shell theory of carbon nanotubes.
	Europ. J. Mech. A/Solids, 42:137–157, 2013.
	
	Brenner_2002
	D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, and S.B.
	Sinnott.
	A second-generation reactive empirical bond order (REBO)
	potential energy expression for hydrocarbons.
	J. Phys. Cond. Mat., 14(4):783, 2002.
	
	Cadelano_2009
	E. Cadelano, P.L. Palla, S. Giordano, and L. Colombo.
	Nonlinear elasticity of monolayer graphene.
	Phys. Rev. Lett., 102:235502, 2009.
	
	Sfyris_2014
	C. Galiotis D. Sfyris, G.I. Sfyris.
	Curvature dependent surface energy for a free standing monolayer
	graphene: Some closed form solutions of the non-linear theory.
	Int. J. nonl. Mech., 67:186–197, 2014.
	
	Sfyris_2014b
	C. Galiotis D. Sfyris, G.I. Sfyris.
	Curvature dependent surface energy for free standing monolayer
	graphene: Geometrical and material linearization with closed form solutions.
	Int. J. Engineering Science, 85:224–233, 2014.
	
	Davini_2014
	C. Davini.
	Homogenization of a graphene sheet.
	Cont. Mech. Thermod., 26(1):95–113, 2014.
	
	Davini_2017
	C. Davini, A. Favata, and R. Paroni.
	A homogenized continuum theory for graphene bending.
	Forthcoming, 2017.
	
	Deng_2016
	S. Deng and V. Berry.
	Wrinkled, rippled and crumpled graphene: an overview of formation
	mechanism, electronic properties, and applications.
	Materials Today, 19(4):197 – 212, 2016.
	
	Favata_2014
	A. Favata, A. Micheletti, and P. Podio-Guidugli.
	A nonlinear theory of prestressed elastic stick-and-spring
	structures.
	Int. J. Eng. Sci., 80:4–20, 2014.
	
	Favata_2016
	A. Favata, A. Micheletti, P. Podio-Guidugli, and N.M. Pugno.
	Geometry and self-stress of single-wall carbon nanotubes and graphene
	via a discrete model based on a 2nd-generation REBO potential.
	J. Elasticity, 125:1–37, 2016.
	
	Favata_2016b
	A. Favata, A. Micheletti, P. Podio-Guidugli, and N.M. Pugno.
	How graphene flexes and stretches under concomitant bending couples
	and tractions.
	Meccanica, in press, 2016.
	
	Favata_2016a
	A. Favata, A. Micheletti, S. Ryu, and N.M. Pugno.
	An analytical benchmark and a mathematica
	program for MD codes: testing LAMMPS on the 2nd
	generation Brenner potential.
	Comput. Phys. Commun., 207:426–431, 2016.
	
	Ferrari2014
	A.C. Ferrari, F. Bonaccorso, V. Fal'ko, K.S. Novoselov, S. Roche,
	P. Bøggild, S. Borini, F.H.L. Koppens, V. Palermo, N.M. Pugno, J.A.
	Garrido, R. Sordan, A. Bianco, L. Ballerini, M. Prato, E. Lidorikis,
	J. Kivioja, C. Marinelli, T. Ryhänen, A. Morpurgo, J.N. Coleman,
	V. Nicolosi, L. Colombo, A. Fert, M. Garcia-Hernandez, A. Bachtold, G.F.
	Schneider, F. Guinea, C. Dekker, M. Barbone, Z. Sun, C. Galiotis, A.N.
	Grigorenko, G. Konstantatos, A. Kis, M. Katsnelson, L. Vandersypen,
	A. Loiseau, V. Morandi, D. Neumaier, E. Treossi, V. Pellegrini, M. Polini,
	A. Tredicucci, G.M. Williams, B. Hee Hong, J.-H. Ahn, J. Min Kim, H. Zirath,
	B.J. van Wees, H. van der Zant, L. Occhipinti, A. Di Matteo, I.A. Kinloch,
	T. Seyller, E. Quesnel, K. Feng, X.and Teo, N. Rupesinghe, P. Hakonen,
	S. R.T. Neil, Q. Tannock, T. Löfwander, and J. Kinaret.
	Science and technology roadmap for graphene, related two-dimensional
		crystals, and hybrid systems.
	Nanoscale, 7(11):4587–5062, 2015.
	
	Goler_2013
	S. Goler, C. Coletti, V. Tozzini, V. Piazza, T. Mashoff, F. Beltram,
	V. Pellegrini, and S. Heun.
	Influence of graphene curvature on hydrogen adsorption: Toward
	hydrogen storage devices.
	Phys. Chem. C, 117(22):11506–11513, 2013.
	
	Hajgato_2012
	B. Hajgató, S. Güryel, Y. Dauphin, J.-M. Blairon, H.E. Miltner,
	G. Van Lier, F. De Proft, and P. Geerlings.
	Theoretical investigation of the intrinsic mechanical properties of
	single- and double-layer graphene.
	J. Phys. Chem. C, 116(42):22608–22618, 2012.
	
	Hartmann_2013
	M.A. Hartmann, M. Todt, F.G. Rammerstorfer, F.D. Fischer, and O. Paris.
	Elastic properties of graphene obtained by computational mechanical
	tests.
	Europhys. Lett., 103(6):68004, 2013.
	
	Huang2006
	Y. Huang, J. Wu, and K. C. Hwang.
	Thickness of graphene and single-wall carbon nanotubes.
	Phys. Rev. B, 74:245413, 2006.
	
	Kim_2012
	S.M. Kim, E.B. Song, S. Lee, J. Zhu, D.H. Seo, M. Mecklenburg, S. Seo, and K.L.
	Wang.
	Transparent and flexible graphene charge-trap memory.
	ACS Nano, 6(9):7879–7884, 2012.
	
	Koskinen_2010
	P. Koskinen and O.O. Kit.
	Approximate modeling of spherical membranes.
	Phys. Rev. B, 82:235420, 2010.
	
	Kudin_2001
	K.N. Kudin, G.E. Scuseria, and B.I. Yakobson.
	c_2F, BN, and C
	nanoshell elasticity from ab initio computations.
	Phys. Rev. B, 64:235406, 2001.
	
	Lindahl_2012
	N. Lindahl, D. Midtvedt, J. Svensson, O.A. Nerushev, N. Lindvall, A. Isacsson,
	and E.E. B. Campbell.
	Determination of the bending rigidity of graphene via electrostatic
	actuation of buckled membranes.
	Nano Lett., 12(7):3526–3531, 2012.
	
	Liu_2007
	F. Liu, P. Ming, and J. Li.
	Ab initio calculation of ideal strength and phonon
	instability of graphene under tension.
	Phys. Rev. B, 76:064120, 2007.
	
	LuJ1997
	J.P. Lu.
	Elastic properties of carbon nanotubes and nanoropes.
	Phys. Rev. Lett., 79:1297–1300, 1997.
	
	Lu_2009
	Q. Lu, M. Arroyo, and R. Huang.
	Elastic bending modulus of monolayer graphene.
	J. Phys. D, 42(10):102002, 2009.
	
	Luhuang2009
	Q. Lu and R. Huang.
	Nonlinear mechanics of single-atomic-layer graphene sheets.
	Int. J. Appl. Mech., 01(03):443–467, 2009.
	
	Pacheco_2014
	A.A. Pacheco Sanjuan, Z. Wang, H.P. Imani, M. Vanević, and
	S. Barraza-Lopez.
	Graphene's morphology and electronic properties from discrete
	differential geometry.
	Phys. Rev. B, 89:121403, 2014.
	
	SakhaeePour_2009
	A. Sakhaee-Pour.
	Elastic properties of single-layered graphene sheet.
	Sol. St. Comm., 149(1–2):91 – 95, 2009.
	
	Scarpa_2010
	F. Scarpa, S. Adhikari, A.J. Gil, and C. Remillat.
	The bending of single layer graphene sheets: the lattice versus
	continuum approach.
	Nanotech., 21(12):125702, 2010.
	
	Scarpa_2009
	F. Scarpa, S. Adhikari, and A. Srikantha Phani.
	Effective elastic mechanical properties of single layer graphene
	sheets.
	Nanotech., 20(6):065709, 2009.
	
	Shi_2012
	X. Shi, B. Peng, N.M. Pugno, and H. Gao.
	Stretch-induced softening of bending rigidity in graphene.
	Appl. Phys. Let., 100(19), 2012.
	
	Tapaszto_2012
	L. Tapaszto, T. Dumitrica, S.J. Kim, P. Nemes-Incze, C. Hwang, and L.P. Biro.
	Breakdown of continuum mechanics for nanometre-wavelength rippling
		of graphene.
	Nat. Phys., 8(10):739–742, 2012.
	
	Tozzini_2011
	V. Tozzini and V. Pellegrini.
	Reversible hydrogen storage by controlled buckling of graphene
	layers.
	Phys. Chem. C, 115(51):25523–25528, 2011.
	
	Tozzini_2013
	V. Tozzini and V. Pellegrini.
	Prospects for hydrogen storage in graphene.
	Phys. Chem., 15:80–89, 2013.
	
	Wei_2013
	Y. Wei, B. Wang, J. Wu, R. Yang, and M.L. Dunn.
	Bending rigidity and Gaussian bending stiffness of
	single-layered graphene.
	Nano Lett., 13(1):26–30, 2013.
	
	Yakobson_1996
	B. I. Yakobson, C. J. Brabec, and J. Bernholc.
	Nanomechanics of carbon tubes: Instabilities beyond linear response.
	Phys. Rev. Lett., 76:2511–2514, Apr 1996.
	
	Zakharchenko_2009
	K. V. Zakharchenko, M. I. Katsnelson, and A. Fasolino.
	Finite temperature lattice properties of graphene beyond the
	quasiharmonic approximation.
	Phys. Rev. Lett., 102:046808, 2009.
	
	Zhang_2011
	D.-B. Zhang, E. Akatyeva, and T. Dumitric ăă.
	Bending ultrathin graphene at the margins of continuum mechanics.
	Phys. Rev. Lett., 106:255503, 2011.
	
	Zhao_2009
	H. Zhao, K. Min, and N. R. Aluru.
	Size and chirality dependent elastic properties of graphene
	nanoribbons under uniaxial tension.
	Nano Lett., 9(8):3012–3015, 2009.
	




