---
author:
- Cesare Davini$^1$
- Antonino Favata$^2$
- Roberto Paroni$^3$
title: |
     **A new material property of graphene:\
    the bending Poisson coefficient**
---

[ $^1$ Via Parenzo 17, 33100 Udine\
<cesare.davini@uniud.it>\
$^2$ Department of Structural and Geotechnical Engineering\
Sapienza University of Rome, Rome, Italy\
<antonino.favata@uniroma1.it>\
$^3$ DADU\
University of Sassari, Alghero (SS), Italy\
<paroni@uniss.it> ]{}

Abstract {#abstract .unnumbered}
========

The in-plane infinitesimal deformations of graphene are well understood: they can be computed by solving the equilibrium problem for a sheet of isotropic elastic material with suitable stretching stiffness and Poisson coefficient $\nu^{(m)}$. Here, we pose the following question: does the Poisson coefficient $\nu^{(m)}$ affect the response to bending of graphene?  Despite what happens if one adopts classical structural models, it does not.  In this letter we show that a new material property, conceptually and quantitatively different from $\nu^{(m)}$, has to be introduced.  We term this new parameter *bending Poisson coefficient*; we propose for it a physical interpretation in terms of the atomic interactions and produce a quantitative evaluation. 

**PACS**: 61.48.Gh (Structure of graphene), 62.20.dj (Poisson’s ratio), 62.20.dq (Other elastic constants)

Introduction
============

The in-plane energy of an isotropic elastic sheet, occupying a two-dimensional region $\Omega$, is $$\label{sheet}
\Uc^{(m)}=\Ec\int_\Omega\nu^{(m)} (\mbox{tr}\, \Eb)^2+(1-\nu^{(m)})|\Eb|^2,$$ where, for $\ub$ the in-plane infinitesimal displacement, $\Eb=1/2 (\nabla\ub+\nabla\ub^T)$ is the strain, $\Ec$ is the *stretching stiffness*, and $\nu^{(m)}$ is the Poisson coefficient.  Instead, the bending energy of the sheet is $$\label{platen}
\Uc^{(b)}=\frac{1}{2}\,\Dc\int_\Omega(\Delta w)^2-2(1-\nu^{(b)})\det\nabla^2 w,$$ where $w$ is the infinitesimal transversal displacement; $\Dc$ and $\nu^{(b)}$ are two parameters, the first of which is the *bending stiffness*.  Energies and are classically deduced from a three-dimensional model from which it results that $\nu^{(b)}=\nu^{(m)}$.

In *Davini*[@Davini_2014], see also references therein, it has been shown that the energy associated to in-plane infinitesimal deformations of graphene has the same form of with $\Ec$ and $\nu^{(m)}$ given in terms of elastic lattice constants. 

Here, we are mainly concerned with the bending behavior of graphene.  At macroscopic level, graphene can be considered as an elastic sheet able to sustain bending, a feature that has attracted a growing interest in the past few years because of the possible technological applications [@Ferrari2014].  In this letter we show that the bending energy associated to infinitesimal displacements of graphene can be recast as in , but the identity $\nu^{(b)}=\nu^{(m)}$, which holds for classical structural models, does not hold.  Indeed, we show that for a graphene sheet $$\nu^{(b)}\ne\nu^{(m)}.$$ We call *bending Poisson coefficient* the new independent parameter $\nu^{(b)}$ . 

![ A sheet of graphene in the shape of a double curvature, anticlastic surface.[]{data-label="fig:bent_graphene"}](bent_graphene)

A major conclusion is therefore that the material parameters $\Dc$ and $\nu^{(m)}$, which are well known in the literature on graphene, are not enough to fully describe the bending behavior.  Indeed, when dealing with the small-deformations, the research has so far mainly focused on the in-plane behavior or on cylindrical bending of graphene; in this latter case, as is clear from , the bending Poisson coefficient has no role because $\det\nabla^2 w=0$.

When double curvature deformations are considered (see Fig.  \[fig:bent\_graphene\]), it is usually supposed that $\nu^{(b)}=\nu^{(m)}$. In *Zhang et al.*[@Zhang_2011], the validity of the energy - is questioned, basically because of the lack of a specific counterpart at the microscopic level; it is there adopted a $\pi$-orbital axis vector model and tight binding computations based on density functional are performed, and some limits on the validity of the standard continuum theory have been described.  Nevertheless, in the simulations carried out just cylindrical bending has been considered, a circumstance that cannot highlight the need of considering a bending Poisson coefficient as already noticed. 

The bending stiffness $\Dc$ has been lengthy investigated.  An account of the literature on the subject would be necessarily incomplete, and the reader is addressed to the very recent review by *Akinwande et al.*[@Akiwande_2017].  In particular, we refer to evaluations based on molecular dynamics computations[@Huang2006; @Lu_2009; @Favata_2016b].  Other techniques, such as density functional theory, tight binding, quantum mechanics-based methods, can be found in the literature [@Lindahl_2012; @Tapaszto_2012; @Kim_2012; @Shi_2012; @Hartmann_2013; @Hajgato_2012; @Wei_2013; @Pacheco_2014].  Direct measurement of $\Dc$ has been challenging for monolayer graphene as well as for other 2D materials.  The value often quoted for the bending modulus of monolayer graphene is $\simeq 1.2$ eV, estimated from the phonon spectrum of bulk graphite [@Nicklow_1972]. * Lindahl et al.*[@Lindahl_2012] determined the bending stiffness of double-layer graphene based on measurements of the critical voltage for snap-through of pre-buckled graphene membranes; the same method was applied to monolayer graphene membranes, yielding a rough estimate with higher uncertainties due to rather limited data point.  The membranal Poisson coefficient $\nu^{(m)}$, on the other hand, has been determined by means of density functional theory[@Kudin_2001; @Wei_2009] and by molecular dynamics simulations[@Brenner_1990; @Brenner_2002; @Stuart2000; @Lindsay_2010]. 

For the new material constant $\nu^{(b)}$, we furnish a nano-scale interpretation and propose an analytical formula induced from the 2nd-generation Brenner potential[@Brenner_2002], which is one of the most used in molecular dynamics simulations of graphene.  This parameter rules the attitude of graphene to form surfaces with double curvature, and is fundamental in designing whatever flexible graphene-based device. 

A MD-induced description of graphene bending {#model}
============================================

Molecular Dynamics (MD) simulations are often employed to study the mechanical behavior of graphene and the 2nd-generation REBO (Reactive Empirical Bond Order) potential is one of the most used empirical potentials.  As is known, this potential was originally developed for hydrocarbons by [*Brenner et al.*]{}[@Brenner_2002], and is able to accommodate up to third-nearest-neighbor interactions. 

In an easy-to-visualize mechanical picture of graphene, the kinematic variables associated with the interatomic bonds involve first, second and third nearest neighbors of any given atom.  In particular, the kinematical variables to be considered are *bond lengths* $l$, *wedge angles* $\vartheta$ and *dihedral angles*.  From *Brenner et al.*[@Brenner_2002], dihedral angles are of two kinds that we here term C and Z and denote by $\operatorname{\overset{(c)}{\Theta}}$ and $\operatorname{\overset{(z)}{\Theta}}$ (see Fig.  \[fig:bonds\]).  The first two contributions are related to the strong covalent $\sigma$-bonds between atoms in one and the same given plane, while the role of dihedral angles is to account for the local coordination of, and the bond angles relative to, two atoms, and are related to the $\pi$-bonds perpendicular to the plane of $\sigma$-bonds. 

![ Kinematic variables: distance $l$, angle $\vartheta$, Z-dihedral angle $\operatorname{\overset{(z)}{\Theta}}$ and C-dihedral angle $\operatorname{\overset{(c)}{\Theta}}$.[]{data-label="fig:bonds"}](bonds)

We consider a harmonic approximation of the stored energy and assume that it is given by the sum of the following terms: $$\label{eq:ENER 1bis}
\begin{array}{l}
\displaystyle \mathcal{U}_\ell^l = \frac{1}{2} \, \sum_{\mathcal {E}} k^l \, (\delta l)^2 ,\\
\displaystyle\mathcal{U}_\ell^\vartheta = \tau_{0} \sum_{\mathcal {W}}\delta \vartheta +  \frac{1}{2} \, \sum_{\mathcal {W}} k^\vartheta \, (\delta\vartheta)^2,\\
\displaystyle\mathcal{U}_\ell^\Theta = \frac{1}{2} \, \sum_{\mathcal {Z}} k^{\Theta}\, (\delta\operatorname{\overset{(z)}{\Theta}})^2+\frac{1}{2} \, \sum_{\mathcal {C}} k^{\Theta}\, (\delta\operatorname{\overset{(c)}{\Theta}})^2. 
\end{array}$$ $\mathcal{U}_\ell^l$, $\mathcal{U}_\ell^\vartheta $ and $\mathcal{U}_\ell^\Theta$ are the energies of the edge bonds, the wedge bonds and the dihedral bonds, respectively; $\tau_0$ is the wedge self-stress, responsible of the cohesive energy.  The role of the wedge self-stress is crucial and its presence is suggested by the shape of the 2nd generation Brenner potential (see *Favata et al.*[@Favata_2016] for details).  The sums in extend to all edges, $\mathcal{E}$, all wedges, $\mathcal{W}$, all Z-dihedra, $\mathcal{Z}$, and all C-dihedra, $\Cc$. The bond constants $k^l$, $k^\vartheta$, and $k^\Theta$ are deduced by making use of the 2nd-generation Brenner potential. 

Here, $\ell$ is the natural length of the C–C bond in the reference configuration; $L^\ell_1$ and $L^\ell_2$ denote the two Bravais lattices realizing the graphene layer (see Fig.  \[fg:LATTICE\], where the nodes of $L^\ell_1$ are represented by blank circles and those of $L^\ell_2$ by black spots); and the attention is restricted to a bounded piece of graphene, that is, to the set of points $\xb^\ell\in L^\ell_1\cup L^\ell_2$ contained in a bounded open set $\Omega\subset \mathbb{R}^2$, cf.  Fig.  \[fg:LATTICE\]. 

In order to express the energies in terms of the nodal displacements, it is expedient to introduce the unit vectors $\pb_i$ ($i=1,2,3$) that, from each node of $L_2^\ell$, point at the nearest neighbor lattice points. 

\[!hbtp\]

![ The hexagonal lattice.[]{data-label="fg:LATTICE"}](hexaglatticenew)

We also introduce the vectors $$\db_1:=\pb_1-\pb_3,\quad\db_2:=\pb_2-\pb_3,\quad\mbox{and}\quad\db_3:=\db_1-\db_2,$$ where $\{\ell\db_1,\ell\db_2\}$ is a couple of lattice vectors that generates the lattice $L_1^\ell$. We approximate the strain measures associated to a change of configuration, described by a displacement field $\ub$, by means of a Taylor expansion truncated at the lowest order that makes the energy quadratic in $\ub$.

We notice that, due to the presence of the wedge self-stress $\tau_0$, also the second order variation of the wedge angles counts.  Indeed, if $\delta \vartheta_i^{(1)}(\xb^\ell)$ and $\delta \vartheta_i^{(2)}(\xb^\ell)$ denote the first and the second order variation of the wedge angles at node $\xb^\ell$, respectively, computations that we omit for brevity yield that $\sum_{i=1}^3 \delta \vartheta_i^{(1)}(\xb^\ell)=0$. Thus, when truncated at the quadratic term, the wedge energy becomes $$\displaystyle\mathcal{U}_\ell^\vartheta =\tau_{0} \sum_{\xb^\ell}\sum_{i=1}^3\delta \vartheta_i^{(2)}(\xb^\ell)+\frac{1}{2}  k^\vartheta \, \sum_{\xb^\ell} \sum_{i=1}^3 (\delta \vartheta_i^{(1)}(\xb^\ell))^2,$$ where, here and below, the sum on $\xb^\ell$ is over the set of points $(L^\ell_2\cup L^\ell_2)\cap\Omega$ unless differently stated. 

Computations also show that $\delta l_i$ and $\delta \vartheta_i^{(1)}$ depend upon the in-plane components of $\ub$, while $\delta \vartheta_i^{(2)}$, $\delta\operatorname{\overset{(c)}{\Theta}}$ and $\delta\operatorname{\overset{(z)}{\Theta}}$ depend upon the out-of-plane component of $\ub$. This yields a splitting of the energy into *membrane* and *bending* parts $$\mathcal{U}_\ell=\mathcal{U}_\ell^{(m)}+\mathcal{U}_\ell^{(b)},$$ with $\mathcal{U}_\ell^{(m)}$ and $\mathcal{U}_\ell^{(b)}$ defined by $$\mathcal{U}_\ell^{(m)}:= \frac{1}{2} k^l \, \sum_{\xb_\ell} \sum_{i=1}^3 (\delta l_i(\xb_\ell))^2 + \frac{1}{2}  k^\vartheta \, \sum_{\xb_\ell} \sum_{i=1}^3 (\delta \vartheta_i^{(1)}(\xb^\ell))^2$$ and $$\begin{aligned}
\mathcal{U}_\ell^{(b)}:=&\tau_{0} \sum_{\xb_\ell}\sum_{i=1}^3\delta \vartheta_i^{(2)}(\xb^\ell)\nonumber\\
&\hspace{3mm}+\frac{1}{2} \, \sum_{\mathcal {Z}} k^{\Theta}\, (\delta\operatorname{\overset{(z)}{\Theta}})^2+\frac{1}{2} \, \sum_{\mathcal {C}} k^{\Theta}\, (\delta\operatorname{\overset{(c)}{\Theta}})^2.\label{bendergy}\end{aligned}$$

From a detailed study of kinematics, it turns out that the expressions giving the in-plane strains $\delta l_i$ and $\delta \vartheta_i^{(1)}$ coincide with those worked out by *Davini*[@Davini_2014] in dealing with the plane deformations of a graphene sheet.  Therefore, that analysis applies and provides a characterization of the continuum limit, for $\ell\to 0$, of the membrane energy $\mathcal{U}_\ell^{(m)}$. So, hereafter we concentrate on the bending energy only. 

The first term on the right hand side of $$\label{selfen1}
\mathcal{U}^{(s)}_\ell:= \tau_0\sum_{\xb^\ell}  \sum_{i = 1}^3\delta \vartheta_i^{(2)} (\xb^\ell),$$ is the *cohesive energy*. 

Fig.  \[fig:cell\_text\] illustrates the meaning of the dihedral angles $\operatorname{\overset{(z)}{\Theta}}_{\pb_i\pb_{i+2}}$, $\operatorname{\overset{(z)}{\Theta}}_{\pb_i\pb_{i+1}}$, $\operatorname{\overset{(c)}{\Theta}}_{\pb_i^+}$, $\operatorname{\overset{(c)}{\Theta}}_{\pb_i^-}$, defined as the angles between the planes spanned by two adjacent segments in colored chain; with these, the other two terms in the energy can be written as $$\begin{aligned}
\label{enZ}
\mathcal{U}^{\mathcal{Z}}_\ell=\frac{1}{2}\, k^{\Theta} \sum_{\xb^\ell\in L^\ell_2\cap\Omega}
\sum_{i=1}^3 &\Bigg(\delta\operatorname{\overset{(z)}{\Theta}}_{\pb_i\pb_{i+2}}(\xb^\ell)\Bigg)^2+\\
&+ \Bigg(\delta\operatorname{\overset{(z)}{\Theta}}_{\pb_i\pb_{i+1}}(\xb^\ell)\Bigg)^2,\nonumber\end{aligned}$$ and $$\label{enC}
\mathcal{U}^{\mathcal{C}}_\ell=\frac{1}{2} k^{\Theta} \, \sum_{\xb^\ell\in L^\ell_2\cap\Omega}
\sum_{i=1}^3 \Bigg(\delta\operatorname{\overset{(c)}{\Theta}}_{\pb_i^+}(\xb^\ell)\Bigg)^2+ \Bigg(\delta\operatorname{\overset{(c)}{\Theta}}_{\pb_i^-}(\xb^\ell)\Bigg)^2,$$ where $\mathcal{U}^{\mathcal{Z}}_\ell$ and $\mathcal{U}^{\mathcal{C}}_\ell$ respectively are the [*Z-dihedral energy*]{} and the [*C-dihedral energy*]{}.  Here, $i, i+1, i+2$ take values in $\{1,2,3\}$ and the sum is to be interpreted modulo 3: for instance, if $i=3$, then $i+1=1$ and $i+2=2$.

![ Left: C-dihedral angles $\operatorname{\overset{(c)}{\Theta}}_{\pb_1^+}$ (green) and Z-dihedral angle $\operatorname{\overset{(z)}{\Theta}}_{\pb_1\pb_2}$ (blue).  Right: C-dihedral angles $\operatorname{\overset{(c)}{\Theta}}_{\pb_1^-}$ (green) and Z-dihedral angle $\operatorname{\overset{(z)}{\Theta}}_{\pb_1\pb_3}$ (blue).[]{data-label="fig:cell_text"}](cell_text)

In and we have taken the elastic constants of both types of dihedral energy to be equal, as implied by the potential adopted.  We notice also that the term $\sum_{i=1}^3\delta \vartheta_i^{(2)} (\xb^\ell)$ in the cohesive energy, cf. ,  is always non-positive, as is made clear in the Appendix.  Hereafter, we assume that $$\tau_0<0\quad\mbox{and}\quad k^{\Theta} >0,$$ so that all the components of the bending energy are non-negative. 

The continuum bending energy
============================

The discrete bending energy $\mathcal{U}_\ell^{(b)}$ depends on out-of-plane displacements defined over the set of points $(L^\ell_1\cup L^\ell_2)\cap\Omega$. We assume that it can be approximated by a continuous energy $\Uc_0^{(b)}$ depending on functions defined over the domain $\Omega$, to be determined by letting the lattice size $\ell$ go to zero so that $(L^\ell_1\cup L^\ell_2)\cap\Omega$ invades $\Omega$.

Lengthy computations, that we omit for the sake of space, yield: $$\begin{aligned}
\delta\operatorname{\overset{(z)}{\Theta}}_{\pb_i\pb_{i+2}}(\xb^\ell)&=\frac{2\sqrt{3}}{3\ell}[w(\xb^\ell+\ell\pb_i-\ell \pb_{i+2})-w(\xb^\ell)\\
&-w(\xb^\ell+\ell \pb_{i})+w(\xb^\ell+\ell \pb_{i+2})]\nonumber\\
&=\frac{2\sqrt{3}}{3\ell}[w(\xb^\ell+\ell\ab_i)-w(\xb^\ell)\\
&-w(\xb^\ell+\ell \ab_{i}+\ell \pb_{i+2})+w(\xb^\ell+\ell \pb_{i+2})]\nonumber\\
&=\frac{2\sqrt{3}}{3\ell} [- \partial^2_{\ab_i\pb_{i+2}}w(\xb^\ell) \ell^2 |\ab_i| |\pb_{i+2}|+o(\ell^2)]\nonumber\\
&=-2 \ell \partial^2_{\ab_i\pb_{i+2}}w(\xb^\ell) +o(\ell),\label{angZ1}\end{aligned}$$ where $w$ is the out-of-plane component of displacement and $\ab_i:=\pb_i-\pb_{i+2}$, $|\ab_i |=\sqrt{3}$; $\partial^2_{\ab_i\pb_{i+2}}$ denotes the second partial derivative in the directions of $\ab_i$ and $\pb_{i+2}$; and $o(\cdot)$ stands for an infinitesimal quantity of order greater than its argument.  With in mind, we rewrite this identity as $$\begin{aligned}
\left(\delta\operatorname{\overset{(z)}{\Theta}}_{\pb_i\pb_{i+2}}(\xb^\ell)\right)^2
=\\
\frac{8\sqrt{3}}9 &\int_{E^\ell(\xb^\ell)}(\partial^2_{\ab_i\pb_{i+2}}w(\xb^\ell))^2\,d\xb+o(\ell^2),\label{angZ2}
\end{aligned}$$ where $E^\ell(\xb^\ell)$ is the hexagon of side $\ell$ centred at $\xb^\ell$, whose area is $\ell^23\sqrt{3}/2$ (see Fig.  \[fg:LATTICE\]).  Similarly, we have: $$\begin{aligned}
\left(\delta\operatorname{\overset{(z)}{\Theta}}_{\pb_i\pb_{i+1}}(\xb^\ell)\right)^2
=\\
\frac{8\sqrt{3}}9 &\int_{E^\ell(\xb^\ell)}(\partial^2_{\bb_i\pb_{i+1}}w(\xb^\ell))^2\,d\xb +o(\ell^2),
\end{aligned}$$ where $\bb_i:=\pb_i-\pb_{i+1}$. By writing $\ab_i$ and $\bb_i$ explicitly, using , and observing that the number of items in the sums has order $1/\ell^2$, for $\ell$ going to zero we find that $$\begin{aligned}
\mathcal{U}^{\mathcal{Z}}_0(w)&:=\lim_{\ell\to 0}\mathcal{U}^{\mathcal{Z}}_\ell=\frac{1}{2}\frac{8\sqrt{3}}9 k^{\Theta}\int_\Omega (\partial^2_{\db_1\pb_{3}}w)^2+ (\partial^2_{\db_1\pb_{1}}w)^2\label{limenZ}\\
&\hspace{-5mm}+(\partial^2_{\db_2\pb_{2}}w)^2
+(\partial^2_{\db_2\pb_{3}}w)^2
+(\partial^2_{\db_3\pb_{1}}w)^2+(\partial^2_{\db_3\pb_{2}}w)^2\,d\xb.\nonumber\end{aligned}$$ Hence, by rewriting the directional derivatives in terms of the partial derivatives with respect to an orthogonal system of coordinates $\{x_1, x_2\}$ we deduce that $$\begin{aligned}
\mathcal{U}^{\mathcal{Z}}_0(w)
&=\frac{1}{2}\frac{5\sqrt{3}}3 k^{\Theta}\int_\Omega  (\Delta w)^2- \frac{8}5 \det \nabla^2 w\,d\xb. 
\label{enZ0}\end{aligned}$$ Analogous expressions can be derived for $\mathcal{U}^{\mathcal{C}}_0(w)$ and $\mathcal{U}^{(s)}_0(w)$, cf.  Appendix.  Then, by summing all contributions up, we find that the continuum limit of the bending energy is given by: $$\begin{aligned}
\label{entot}
\mathcal{U}^{(b)}_0(w)
&=\frac 12 \Bigg(\frac{7\sqrt{3}}3 k^{\Theta}-\frac{\tau_0}2\Bigg) \int_\Omega( \Delta w)^2\nonumber\\
&\hspace{0.5cm}-2\,\frac{16 k^\Theta}{14 k^\Theta-\sqrt{3}\tau_0}\det \nabla^2 w\,d\xb\nonumber\\
&=\frac 12 \Dc  \int_\Omega ( \Delta w)^2-2(1-\nu^{(b)})\det \nabla^2 w \,d\xb,$$ where we have set $$\label{nub}
\nu^{(b)}=1-\frac{16 k^\Theta}{14 k^\Theta-\sqrt{3}\tau_0}.$$ We call $\nu^{(b)}$ the *bending Poisson coefficient*. 

In sheets of conventional materials, the coefficient $\nu^{(b)}$ appearing in eq.  is the membranal Poisson coefficient $\nu^{(m)}$. For graphene this is no longer so.  In fact from the continuum in-plane energy of graphene one finds: $$\label{numm}
\nu^{(m)}=\frac{k^l\ell^2-6k^\vartheta}{k^l\ell^2+18k^\vartheta}\neq\nu^{(b)},$$ (see the paper by *Davini*[@Davini_2014] for a detailed computation).  Besides the quantitative evaluations, detailed in the next section, it is worth noticing that the two coefficients are conceptually different, as the nano-scale interpretation given by and clearly reveals.  In fact, the membranal Poisson coefficient depends upon

(i) the reluctance of two arbitrary atoms to change their distance, captured by $k^l$;

(ii) the reluctance of three atoms to change their mutual angle, captured by $k^\vartheta$.

The bending Poisson coefficient, on the contrary, is determined by

(i) the reluctance of four arbitrary atoms to change the dihedral angles they form, captured by $k^\Theta$,

(ii) the amount of the cohesive energy, related to $\tau_0$.

Equation also confirms that the origin of the bending stiffness $$\label{bendstiff}
\Dc=\frac{7\sqrt{3}}3 k^{\Theta}-\frac{\tau_0}2$$ is twofold: a part depends on the dihedral contribution, and a part on the cohesive energy (cf. * Favata et al.*[@Favata_2016b]).  Computations following from the quantitative evaluation of $k^{\Theta}$ and $\tau_0$ in the next section show that the two terms contribute to $\mathcal D$ almost equally. 

Results
=======

In order to obtain proper quantitative evaluations of the bending Poisson coefficient, let us recall the 2-nd generation Brenner potential[@Brenner_2002].  The binding energy $V$ of an atomic aggregate is given as a sum over nearest neighbors: $$\label{V}
V=\sum_i\sum_{j<i} V_{ij}\,;$$ the interatomic potential $V_{ij}$ is given by $$\label{Vij}
V_{ij}=V_R(r_{ij})+b_{ij}V_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}$, which depends on the angles $\vartheta_{ijk}$ and the dihedral angles $\Theta_{ijkl}$. The values of the constant $k^\Theta$ can be then deduced by deriving twice the potential, and computing the result in the ground state, where $r_{ij}=\ell$, $\theta_{ijk}=2/3\pi$, $\Theta_{ijkl}=0$. Moreover, from *Favata et al.*[@Favata_2016], we take the value of the selfstress $\tau_0$: $$\label{tau0}
\tau_0=-0.2209\;\textrm{nN nm}=-1.3787\;\textrm{eV}.$$ It turns out that the bending Poisson coefficient is: $$\nu^{(b)}=0.419.$$ With the values of the constants $k^l$ and $k^\vartheta$ deduced from the 2nd-generation Brenner potential, we find that the membranal Poisson coefficient is $$\nu^{(m)}=0.397.$$

This estimate agrees with the values reported in the literature.  For instance, with results obtained by MD simulations that adopt the 2nd-generation Brenner potential [@Arroyo2004; @Akiwande_2017].  For a further check, we can compare the value of the bending stiffness obtained from with others found in the literature.  Indeed, from we get that $\Dc=1.4022$ eV, which is in good agreement with the results given in references[@LuJ1997; @Wei_2013; @Akiwande_2017], although they are deduced by completely different approaches. 

Conclusions
===========

In summary, based on the description of atomic interactions provided by the 2nd-generation Brenner potential, we have found a new material parameter that influences the bending behavior of graphene.  This parameter, called bending Poisson coefficient, is conceptually and quantitatively different from the already known membranal Poisson coefficient.  We have proposed an analytical evaluation in terms of atomistic quantities, revealing its nano-scale physical sources. 

Acknowledgments {#acknowledgments .unnumbered}
===============

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

Appendix: Deduction of $\mathcal{U}^{\mathcal{C}}_0$ and $\mathcal{U}^{(s)}_0$ {#appendix}
==============================================================================

For completeness, here we sketch the deduction of the expressions of $\mathcal{U}^{\mathcal{C}}_0$ and $\mathcal{U}^{(s)}_0$.

By Taylor expansion theorem, we calculate the change of the $\mathcal{C}$-dihedral angle $\delta\operatorname{\overset{(c)}{\Theta}}_{\pb_i^+}(\xb^\ell)$: $$\begin{aligned}
\label{thC+}
    \delta\operatorname{\overset{(c)}{\Theta}}_{\pb_i^+}(\xb^\ell)&=\frac{2\sqrt{3}}{3\ell}[2w(\xb^\ell)-w(\xb^\ell+\ell \pb_{i+1})\nonumber\\
    &\hspace{10mm}+w(\xb^\ell+\ell\pb_i-\ell \pb_{i+2})-2w(\xb^\ell+\ell \pb_{i})]\nonumber\\
    &={2\ell}\partial^2_{\pb_i\pb_i^\perp}w(\xb^\ell)+o(\ell),
    \end{aligned}$$ where $\pb_i^\perp:=\eb_3\times\pb_i$. Similarly, we find that $$\label{thC-}
    \delta\operatorname{\overset{(c)}{\Theta}}_{\pb_i^-}(\xb^\ell)={2\ell}\partial^2_{\pb_i\pb_i^\perp}w(\xb^\ell)+o(\ell).$$

With in mind, we write the identities $$\begin{array}{l}
    \displaystyle\left(\delta\operatorname{\overset{(c)}{\Theta}}_{\pb_i^+}(\xb^\ell)\right)^2=\frac{8\sqrt{3}}{9}\int_{E^\ell(\xb^\ell)}(\partial^2_{\pb_i\pb_i^\perp}w(\xb^\ell))^2\,d\xb +o(\ell^2)\\
    \displaystyle\left(\delta\operatorname{\overset{(c)}{\Theta}}_{\pb_i^-}(\xb^\ell)\right)^2=\frac{8\sqrt{3}}{9}\int_{E^\ell(\xb^\ell)}(\partial^2_{\pb_i\pb_i^\perp}w(\xb^\ell))^2\,d\xb +o(\ell^2)
    \end{array}$$ and deduce that $$\mathcal{U}^{\mathcal{C}}_\ell=\frac{8\sqrt{3}}{9} k^{\Theta} \sum_{\xb^\ell\in L^\ell_2\cap\Omega}\int_{E^\ell(\xb^\ell)} \sum_{1=1}^3(\partial^2_{\pb_i\pb_i^\perp}w(\xb^\ell))^2\,d\xb+o(1).$$ By passing to the limit we find $$\label{limenC2}
    \mathcal{U}^{\mathcal{C}}_0(w):=\lim_{\ell\to 0}\mathcal{U}^{\mathcal{C}}_\ell=\frac{8\sqrt{3}}{9} k^{\mathcal{C}}\int_\Omega \sum_{i=1}^3 \Big(\partial^2_{\pb_i\pb_i^\perp}w\Big)^2\,d\xb.$$

We now compute the limit of the cohesive energy . 

We focus on $\xb^\ell\in L_2^\ell$, since $\xb^\ell\in L_1^\ell$ produces the same result with the same steps.  Computations show that $\delta \vartheta_i^{(2)} (\xb^\ell)$ has the form $$\label{dth2}
    \sum_{i=1}^3\delta\vartheta_i^{(2)}(\xb^\ell)=-\frac{3\sqrt{3}}{\ell^2}\Bigg(\frac 13 \sum_{i=1}^3w(\xb^\ell+\ell\pb_i)-w(\xb^\ell)\Bigg)^2.$$ So, in particular, $\delta \vartheta_i^{(2)} (\xb^\ell)$ is non-positive definite. 

By using Taylor’s expansion in , we find $$\begin{aligned}
\label{sumthi}
    \sum_{i=1}^3\delta\vartheta_i^{(2)}(\xb^\ell)&=-\frac{\sqrt{3}}{3} {\ell^2}\Big( \partial^2_{\pb_1 \pb_1}w(\xb^\ell)\nonumber\\
    &\hspace{5mm}+\partial^2_{\pb_2 \pb_2}w(\xb^\ell)+\partial^2_{\pb_1 \pb_2}w(\xb^\ell)\Big)^2+o(\ell^2),
    \end{aligned}$$ that can be written as $$\begin{aligned}
\label{sumthi2}
    \sum_{i=1}^3\delta\vartheta_i^{(2)}(\xb^\ell)&=-\frac{4}{9}\int_{T^\ell(\xb^\ell)}\Big( \partial^2_{\pb_1 \pb_1}w(\xb^\ell)+\partial^2_{\pb_2 \pb_2}w(\xb^\ell)\nonumber\\
    &\hspace{20mm}+\partial^2_{\pb_1 \pb_2}w(\xb^\ell)\Big)^2\,dx+o(\ell^2),
    \end{aligned}$$ where $T^\ell(\xb^\ell)$ is a triangle centered at $\xb^\ell\in L_1^\ell\cup L_2^\ell$, with vertices in the center of the hexagonal cells having $\xb^\ell$ in common, see Fig.  \[fig:triangle\].  Note that the area of $T^\ell(\xb^\ell)$ is $\frac{3\sqrt{3}}4\ell^2$.

![ Triangulation $T^\ell(\xb^\ell)$[]{data-label="fig:triangle"}](triangle)

With this expression the cohesive energy takes the form $$\begin{aligned}
    \mathcal{U}^{(s)}_\ell
    &=-\frac{4}{9}\tau_0\sum_{\xb^\ell\in (L_1^\ell\cup L_2^\ell)\cap\Omega}   
    \int_{T^\ell(\xb^\ell)}\Big( \partial^2_{\pb_1 \pb_1}w(\xb^\ell)\nonumber\\
    &\hspace{10mm}+\partial^2_{\pb_2 \pb_2}w(\xb^\ell)+\partial^2_{\pb_1 \pb_2}w(\xb^\ell)\Big)^2\,dx+o(1).  \nonumber 
    \end{aligned}$$ Thence, by passing to the limit, we get $$\begin{aligned}
\label{limentheta2}
    \mathcal{U}^{(s)}_0(w):=&\lim_{\ell\to 0}\mathcal{U}^{(s)}_\ell\nonumber\\
    =&-\frac{4}9\tau_0 \int_\Omega ( \partial^2_{\pb_1 \pb_1}w+\partial^2_{\pb_2 \pb_2}w+\partial^2_{\pb_1 \pb_2}w\Big)^2\, d\xb. 
    \end{aligned}$$

[10]{}

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.u.  A review on mechanics and mechanical properties of 2d materials—graphene and beyond. ,  13:42 – 77, 2017. 

M. Arroyo and T. Belytschko.  Finite crystal elasticity of carbon nanotubes based on the exponential cauchy-born rule. ,  69:115415, 2004. 

D.W. Brenner.  Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films. ,  42(15), 1990. 

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 () potential energy expression for hydrocarbons. ,  14(4):783, 2002. 

C. Davini.  Homogenization of a graphene sheet. ,  26(1):95–113, 2014. 

A.C. Ferrari et al. . ,  7(11):4587–5062, 2015. 

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 potential. ,  125:1–37, 2016. 

A. Favata, A. Micheletti, P. Podio-Guidugli, and N.M. Pugno.  How graphene flexes and stretches under concomitant bending couples and tractions. ,  in press, 2016. 

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. ,  116(42):22608–22618, 2012. 

M.A. Hartmann, M. Todt, F.G. Rammerstorfer, F.D. Fischer, and O. Paris.  Elastic properties of graphene obtained by computational mechanical tests. ,  103(6):68004, 2013. 

Y. Huang, J. Wu, and K. C. Hwang.  Thickness of graphene and single-wall carbon nanotubes. ,  74:245413, 2006. 

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. ,  6(9):7879–7884, 2012. 

K.N. Kudin, G.E. Scuseria, and B.I. Yakobson.  ${\uppercase{c}}_{2}\uppercase{F},$ , and nanoshell elasticity from *ab initio* computations. ,  64:235406, 2001. 

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. ,  12(7):3526–3531, 2012. 

L. Lindsay and D. A. Broido.  Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene. ,  81:205441, May 2010. 

J.P. Lu.  Elastic properties of carbon nanotubes and nanoropes. ,  79:1297–1300, 1997. 

Q. Lu, M. Arroyo, and R. Huang.  Elastic bending modulus of monolayer graphene. ,  42(10):102002, 2009. 

R. Nicklow, N. Wakabayashi, and H. G. Smith.  Lattice dynamics of pyrolytic graphite. ,  5:4951–4962, Jun 1972. 

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. ,  89:121403, 2014. 

X. Shi, B. Peng, N.M. Pugno, and H. Gao.  Stretch-induced softening of bending rigidity in graphene. ,  100(19), 2012. 

S. J. Stuart, A. B. Tutein, and J. A. Harrison.  A reactive potential for hydrocarbons with intermolecular interactions. ,  112:6472, 2000. 

L. Tapaszto, T. Dumitrica, S.J. Kim, P. Nemes-Incze, C. Hwang, and L.P. Biro. . ,  8(10):739–742, 2012. 

X. Wei, B. Fragneaud, C.A. Marianetti, and J.W. Kysar.  Nonlinear elastic behavior of graphene: Ab initio calculations to continuum description. ,  80:205407, Nov 2009. 

Y. Wei, B. Wang, J. Wu, R. Yang, and M.L. Dunn.  Bending rigidity and gaussian bending stiffness of single-layered graphene. ,  13(1):26–30, 2013. 

D.-B. Zhang, E. Akatyeva, and T. Dumitric ăă.  Bending ultrathin graphene at the margins of continuum mechanics. ,  106:255503, Jun 2011.