Estimating the effective bending rigidity of multi-layer graphene

We present a novel analytical prediction for the effective bending rigidity γ eff of multi–layer graphene sheets. Our approach involves using a variational model to determine the folding conformation of multi–layer graphene sheets where the curvature of each graphene layer is taken into account. The Lennard–Jones potential is used to determine the van der Waals interaction energy per unit area and the spacing distance between graphene layers. The mid–line of the folded multi–layer graphene is described by a solution derived in previous work for folded single– and multi–layer graphene. Several curves are obtained for the single–layer solution using different values of the bending rigidity γ, and compared to the mid–line of the folded multi–layer graphene. The total area between these curves and the mid–line is calculated, and the value of γ eff is determined by the single–layer curve for which this area is minimized. While there is some disagreement in the literature regarding the relationship between the bending rigidity and the number of layers, our analysis reveals that the bending rigidity of multi–layer graphene follows an approximate square–power relationship with the number of layers N, where N < 7. This trend is in line with theoretical and experimental studies reported in the literature.


Introduction
Graphene has been the subject of many theoretical and experimental studies due to its unique properties [1][2][3] and potential applications [4][5][6]. It may be considered as the main building material for all other carbon nanostructures since many carbon allotropes may be considered to be constructed from the planar graphene sheet [7]. For example, carbon nanotubes can be thought of as rolled up sheets of graphene [8], and fullerenes can be considered to be wrapped up sheets of graphene [9]. Additionally, novel carbon nanostructures may be designed by reconfiguring graphene in ways that also achieve new properties of these structures [10,11]. Thus we may conclude that the bending rigidity of graphene is often a key factor in determining the properties of novel carbon nanostructures. The bending rigidity is also shown to directly impact the properties of graphene [12] as well as various morphologies of graphene, such as in graphene folds and ripples [13,14].
Previous studies have estimated the bending rigidity of single-layer graphene with empirical potential methods [15,16] and ab initio calculations [17,18]. These studies suggest that the bending rigidity of single-layer graphene ranges from 0.83 eV to 1.61 eV, and they are analysed in detail by Wei et al [19]. The experimental study, performed by Nicklow et al [20], found the bending rigidity of single-layer graphene to be 1.2 eV which fits the reported range very well. While the bending rigidity of single-layer graphene is relatively well characterised, there less certainty regarding the bending rigidity of multi-layer graphene.
Although the bending rigidity of multi-layer graphene has not been dealt with in as much depth as singlelayer graphene, some theoretical approaches have been proposed. Examining self-folding conformations on a flat substrate, Chen et al [21] report that the bending rigidity of 2-to 6-layer graphene follows a square-power relationship with its thickness. Using continuum theory and atomistic simulations, Shen and Wu [22] show that the bending rigidity is linearly proportional to the number of layers when the number of layers exceeds five. The fold of unsupported multi-layer graphene has been modelled previously by Meng et al [23], where they utilized their earlier solution of folded single-layer graphene. In that study, the bending rigidity of multi-layer graphene is modelled by multiplying the bending rigidity of the single-layer graphene by the number of layers, and this represents a linear relationship between the bending rigidity and the number of layers. In the present paper we use a mathematical model for folding unsupported multi-layer graphene based on the calculus of variations, and then compare it to a solution derived earlier for folded single-layer graphene [13]. By changing the bending rigidity of the single-layer graphene, we determine an effective bending rigidity γ eff as a function of the number of layers N.
The remainder of this paper is divided into five sections. In the next section, we present a general mathematical formulation to model the folding conformation of graphene. In section 3 we derive parametric solutions for folded multi-layer graphene taking into account the curvature of each graphene layer. Following this, in section 4, we provide simplified solutions for the mid-line of the fold by considering only the curvature of the mid-line while treating the bending rigidity as a parameter somewhat like a fitting parameter. In section 5 we investigate the effective bending rigidity of multi-layer graphene and compare our results with earlier findings reported in the literature. A summary and some concluding remarks are presented in section 6. The van der Waals (vdW) interaction energy is calculated based on the Lennard-Jones (LJ) potential, and detailed derivations for this parameter are included in appendix.

Model formulation
Recently, variational calculus was employed by the present authors to formulate a continuous model for the folding conformation of multi-layer graphene sheets supported on a substrate [24]. In the present paper, similar techniques to those of [24] are adopted to model the folding conformation of unsupported multi-layer graphene sheets and investigate the effective bending rigidity γ eff . For example, a representative plot for the geometry of folded 3-layer graphene sheets is shown in figure 1. The model is constructed based upon the mid-line of a multi-layer graphene stack which is the middle graphene layer when the number of layers is odd (N = 2n + 1), or a notional layer when the number of layers is even (N = 2n), for some integer n. Assuming a translational symmetry in the fold direction, a two-dimensional problem is considered in this work. Also the fold is assumed to possess a reflective symmetry in the x-axis, and therefore we need only analyse the upper half of the mid-line which is the solid curve in figure 1.
Furthermore, the upper half of the mid-line is divided into three curves according to the sign of the line curvature κ m and line gradient dy/dx. The first curve is bounded by the points (x 0 , y 0 ) and (x 1 , y 1 ) with a strictly negative curvature and is denoted by C 0 . The second curve is bounded by the points (x 1 , y 1 ) and (x 2 , y 2 ) with a strictly positive curvature and is denoted by C 1 . The flat region is denoted by C 2 , from the point (x 2 , y 2 ) to the point (x 3 , y 3 ) and has zero curvature. Also we use C to denote the concatenation of all these three curves.
The folding conformation of graphene is determined based on the elastic bending energy of the graphene and the vdW interaction energy between graphene layers [13]. The elastic bending energy is assumed to be proportional to the square of the curvature. We note that the curvature vanishes in C 2 , and hence the elastic bending energy only originates in the folded regions of the curve C. Thus, we express the elastic energy E e mathematically as where γ is a positive constant denotes the bending rigidity of graphene, s is the arc length, and κ 2 is the total squared curvature. Additionally, we approximate the contribution of the vdW interaction energy by assuming that it dominates only in the flat region that is the curve C 2 , and hence it is modelled as where ò is a positive constant denotes the vdW interaction energy between graphene layers. Hence, the total energy for this problem is given by Furthermore, we assume that the graphene has a fixed length, and therefore an isoperimetric constraint is imposed on the total arc length of the curve C, that is These considerations are identical to those considered in our recent work [24], where the total energy functional to be minimized is derived and given by where s 2 is the arc length at the point (x 2 , y 2 ). As shown in figure 1, if we assume that the fold of the mid-line starts at the origin of the plane, then the functional (3) is subject to the following boundary conditions and the value of x 2 is to be determined from the natural boundary condition. It is anticipated that as we increase the number of layers N, the vdW interaction energy ò increases, but the spacing distance between adjacent layers 2δ g decreases. Therefore, the LJ potential is employed to calculate values for δ g and ò for different number of layers N. Table 1 presents the values of δ g and ò used in this work, which are obtained from the LJ potential. Details for the derivation of these parameters are included in appendix. In the next two sections we prescribe two different quantities for the the total squared curvature κ 2 based on different assumptions, and for each case we present the corresponding parametric solution for the mid-line.

Multi-layer solution
In this section we propose to model the folding conformation of N graphene layers stacked on top of each other. As reported by Wei et al [19], the bending rigidity of single-layer graphene is found to be in the range of 0.83-1.61 eV, and for the purpose of the analysis here, we prescribe γ = 1.0 eV for the bending rigidity of each graphene layer. Also we account for the curvature of each graphene layer and calculate the total squared curvature κ 2 as and where κ m denotes the mid-line curvature and is given by where dots here denote differentiation with respect to s. Moreover, we assume that the layers maintain an equal distance 2δ g from one another. Thus, if the radius of curvature of the mid-line is given by then the radius of curvature of the j th outer/inner curve is given by when N is odd, or ν = 2j − 1 when N is even. Consequently, the corresponding j th outer/inner line curvature is given by We now use the definition of the outer/inner line curvature given by (6) to simplify (5). For each pair of inner and outer curves a distance νδ g away from the mid-line, we may combine their squared curvatures to simplify (5) as follows As shown in [24], a Taylor series expansion may be used to derive three different approximations of the total squared curvature κ 2 , and here we consider the three-term approximation, which is given by Substituting (7) into (3) gives the functional to be minimized subject to the boundary conditions (4).
Variational calculus and nondimensionalization methods are then employed in [24] to minimize this functional. Also, a change of variables is made such that where S is the nondimensional arc length and primes denote derivatives with respect to θ. Then the obtained expression of the nondimensionalized curvature k q m ( ) is given by , and 2 5 5 where β Y is an arbitrary constant arising when integrating the Euler-Lagrange equation.
Taking into account the boundary conditions, the mid-line is represented by a parametric solution arrived at by integrating the following pair of differential equations q k q sin . or the curve C 1 , we seek a continuous solution at the point (X 1 , Y 1 ) where k m changes sign, so that to satisfy this constraint we change the signs of the parametric solutions which yields that the curve C 1 is given by sin d , Parametric solutions for the j th outer/inner curve may be determined using where X(θ) and Y(θ) represent the parametric solution of the mid-line, and ν = 2j when N is odd, or ν = 2j − 1 when N is even, for Î ¼ j n 1, , { }. As in [24], multiplying these solutions by the scaling factor a g =  N m recovers the dimensional solutions for the folded curve of C. The unknown parameter β Y is determined using the boundary condition y 2 = Nδ g , and hence our solution for the mid-line of the fold is fully determined.

Single-layer solution
In this section we develop a simplified model that can be used to obtain the conformation of the mid-line of the fold. Here we only consider the curvature of the mid-line, and thus the total squared curvature κ 2 is reduced to (3) gives the functional to be minimized subject to the boundary conditions (4). We note that a solution to this style of problem has been previously derived by Cox et al [13] to model the fold of a singlelayer graphene sheet. Thus we exploit that solution to generate the folding conformation of the mid-line, although here we treat the parameter γ as a fitting parameter rather than a material property. The substitution θ = 2f − π/2 is used where the obtained solution is given in terms of two functions, namely where F(f, k) and E(f, k) represent the incomplete elliptic integrals of the first and second kind, respectively, with the elliptic , β is an arbitrary constant of integration, and α s = ò/γ as defined in [13]. For C 0 , the curve with negative curvature, the parametric solution is given by For C 1 , the curve with positive curvature, the parametric solution is given by The parameter β may be determined by solving numerically, and hence fully determine our solution for the mid-line.

Results
In this section, we investigate the effective bending rigidity γ eff of multi-layer graphene by examining the folding conformation of the mid-line of the graphene stack. Figure 2 illustrates our approach of reporting the effective bending rigidity γ eff of 2-layer graphene. As in figure 2(a), the predicted folding profile of the mid-line, obtained from the solutions presented in section 3, is compared to several curves, obtained from the solution presented in section 4, with a range of values for the parameter γ. The area between each curve and the mid-line (A) is calculated as shown in figure 2(b), and the value of γ eff is prescribed such that this area is minimized ( = A A min ). Representative plots for the predicted folding profiles of multi-layer graphene with γ = 1.0 eV, and the predicted folding profile of single-layer graphene with γ = γ eff are shown in figure 2(c) and figure 3. Also, detailed information for our prediction of the effective bending rigidity γ eff of multi-layer graphene are included in table 2. We comment that these values of γ eff are reported according to the folding conformations obtained form the multi-layer solution where we use γ = 1.0 eV for the bending rigidity of each layer in the graphene stack. Also, we notice that the effective bending rigidity for multi-layer graphene follows an approximate square-power relationship with the number of layers N, and hence we may express the effective bending rigidity in terms of γ and N as where this relationship is also confirmed by a log-log plot as shown in figure 4(a). It is evident from figure 4(b) that the results obtained from this model for the effective bending rigidity of 2-to 6-layer graphene are slightly higher than the existing results reported by Chen et al [21]. However, our results confirm and capture their idea of the approximate square-power relationship very well. Furthermore, the obtained values of γ eff for 2-layer, 3layer, 5-layer, and 7-layer graphene are in line with the experimental results reported by Han et al [25] which are 2.6-5.8 eV, 3.7-12.2 eV, 26.0 eV, and 12.0-53.0 eV, respectively. In contrast to the assumed linear relationship [23], this model shows that simply using the single-layer solution, with the bending rigidity being multiplied by the number of layers, would not lead to an accurate folding conformation of multi-layer graphene. We also mention that these obtained values of γ eff are based on the choice of γ = 1.0 eV for the multi-layer solution, and any lower (or higher) value of γ, within the physically meaningful range of 0.83-1.61 eV [19], would lead to better (or worse) agreement with the experimental results   [21], but the scope of this work is to show an analytical approach for estimating the effective bending rigidity of multi-layer graphene. In table 3 we summarize the obtained (or assumed) relationships between the bending rigidity of multi-layer graphene and the number of layers N, as reported in the literature. Finally, we comment that as the number of layers N increases, the minimum area A min also increases, and hence we may conclude that the approach we employ here is only suitable for predicting γ eff when the number of layers N is relatively few (our results suggest that N 6). Our approach becomes increasingly inapplicable for larger N and we comment that other researchers report that the bending rigidity exhibits an approximately linear proportionality to the number of layers when N 6 [22].

Summary
We present a novel analytical approach to predict the effective bending rigidity of multi-layer graphene. Exploiting an earlier work by the present authors, we derive a mathematical model to predict the folding conformation of unsupported multi-layer graphene sheets, where the values for the parameters ò and δ g are calculated based on the LJ potential. Then, the predicted folding profile of multi-layer graphene sheets is compared to several curves obtained from the single-layer solution for a range of values of the bending rigidity. The effective bending rigidity is determined such that the area bounded by the mid-line and a curve obtained from the single-layer solution is minimized. For N < 7, we find that the effective bending rigidity approximately follows a square-power relationship as given by (9) which is confirmed by a log-log plot presented in figure 4(a). This relationship correlates favourably with some earlier work such as that given by Chen et al [21], and the values of the effective bending rigidity obtained here are also consistent with experimental results [25].
As we consider the perfect graphene with uniform thickness distribution in this work, we mention that the geometry variation with additive-layer segments or material impurity with residual contaminants, which are not accounted for in this work, may lead to the imperfect graphene [26,27]. The main difference is to account for the elastic energy of the impure region for the imperfect graphene, but this elastic energy vanishes in the perfect graphene. While the perfect graphene yields a unique folded configuration, there might be different geometries for folded imperfect graphene depending on the length of the impure region [26,27]. This leads to a potential future work by applying the current model to study the folding conformations of the imperfect graphene and investigate the change of the effective bending rigidity.  The vdW interaction energy and the spacing distance between graphene layers are calculated based on the LJ potential. We apply similar considerations to those used in our recent paper [24], where we calculated the vdW interaction energy and the equilibrium distances for multi-layer graphene supported on a substrate. Since we model the fold of unsupported multi-layer graphene, the calculations are simpler here because we only need to account for the interactions between graphene layers. For different number of graphene layers, we approximate the vdW interaction energy ò as where ò 2N denotes the interaction energy between 2N graphene layers and ò N denotes the interaction energy between N graphene layers, respectively. The interaction energy between any pair of graphene layers ò g−g is calculated using a formula derived by Baowan et al [28], given by where η C = 0.3812 Å −2 denotes the surface density of carbon atoms on a sheet of graphene, and A = 4ξσ 6 , and B = 4ξσ 12 are positive constants of the attraction and repulsion, respectively. The parameters ξ = 0.00455 eV and σ = 3.431 Å are the LJ parameters and their values are taken from [29]. It is anticipated that as the number of layers changes, the equilibrium distance δ g changes, and thus it is also considered in this calculation. The value of δ g = δ is determined for different number of layers N independently such that dò N /dδ = 0 is satisfied. The values of the vdW interaction energy ò and the spacing distance δ g obtained from this calculation are presented in table 1.