Multi–layer graphene folds supported on a substrate: a variational model

A mathematical model is developed to study the folding behaviour of multi–layer graphene sheets supported on a substrate. The conformation of the fold is determined from variational considerations based on two energies, namely the graphene elastic energy and the van der Waals (vdW) interaction energy between graphene layers and the substrate. The model is nondimensionalized and variational calculus techniques are then employed to determine the conformation of the fold. The Lennard–Jones potential is used to determine the vdW interaction energy as well as the graphene–substrate and graphene–graphene spacing distances. The folding conformation is investigated under three different approximations of the total line curvature. Our findings show good agreement with experimental measurements of multi–layer graphene folds from the literature.


Introduction
Two-dimensional (2D) carbon allotropes attract widespread interest due to their unique properties which make them key components in many applications such as nanoelectronic devices, and biomedical applications [1,2]. One such material is graphene, which is the most commonly investigated allotrope of carbon with structures including fullerenes, carbon nanotubes, carbon nanohorns, and others. These structures have superior electronic [3], mechanical [4], and thermal properties [5,6]. Particularly, graphene has many applications in the field of biomedicine including transport systems, sensors, tissue engineering and biological agents [2]. Nonplanar structures of graphene widen the possible range of applications that exploit these unique properties. Some techniques, such as applying a mechanical force [7], have been used to fold 2D materials into specific structures in order to achieve desirable properties. For example, the process of constructing graphene origami, a programmable nanoscale building block, involves graphene folding [8][9][10][11].
Although graphene is known in one sense as the strongest material ever discovered [4], it may also be folded and unfolded repeatedly via the tip of an atomic force microscope [12]. Both experimental and theoretical studies show that the conformation of folded graphene is determined by the bending rigidity and the vdW interaction strength in the flat region [13,14]. If the vdW energy in the flat region exceeds the bending energy in the curved region, then the fold would be energetically stable, and vice versa. Therefore the stability of the fold depends on the length of the graphene sheet.
A small-deformation model, based on conjugate gradient methods, is developed by Cranford et al [15] to calculate the bending stiffness of graphene sheets and determine the critical length for which the stable fold occurs. Also, these authors develop a coarse-grained model which confirms their results obtained from the theoretical model. Meng et al [16] employ a finite-deformation beam theory, combined with molecular dynamics simulations, to study the self-folding of single-layer graphene. A comparison between small-and finite-deformation models is made in [16], which shows that the folding profile is accurately predicted using finite-deformation models. An exact theoretical model is developed by Lu et al [17] to investigate the geometrical characteristics of selffolded carbon nanotube, and the shape of the folded edge is given in analytical forms involving elliptic integrals. This theoretical model is largely equivalent to the variational approach used by Cox et al [13] to study the folding behaviour of a single-layer graphene sheet. The folding conformation is determined in [13], through minimizing the potential energy functional using the calculus of variations, and the predicted folding profile perfectly matches the experimental images taken by Lopez-Bezanilla et al [14]. The variational approach is also proven successful in predicting the conformations of folded [18], rippled [19], and wrinkled single-layer graphene [20]. As multi-layer graphene is composed of stacked single-layer sheets of graphene, we propose a similar variational approach to study the folding conformation of multi-layer graphene sheets supported on a substrate taking into account some factors which may affect the folding conformation. In each approximation of the total line curvature, a variational approach is utilized to derive expressions for the curvature, and then parametric solutions are obtained based on these expressions.
Various factors may have a substantial influence on the folding conformation of multi-layer graphene sheets. For example, Lopez-Bezanilla et al [14] outline several differences in the folding profiles depending on the number of graphene layers. Therefore we take this into account by calculating the total curvature for a number of lines each denoting an individual graphene layer. Furthermore, the vdW interactions between neighbouring graphene sheets, as well as between each graphene sheet and the substrate when it is supported, are expected to affect the folding mechanics of multi-layer graphene sheets. The effects of vdW interactions have not previously been incorporated into a mathematical model, partly due to limited experimental studies reported in the literature. In the present work we take into consideration the effects of vdW interactions on the folding conformation of multi-layer graphene, that is supported on a substrate, and compare our findings with experimental data in the literature.
In this paper we propose a mathematical model for the fold of multi-layer graphene located on a substrate. In the following section, we formulate the general mathematical model and approximate the total squared curvatures using a Taylor series expansion. In section 3, we nondimensionalize the model and apply variational calculus methods to derive three expressions for the line curvature each of which is under a different approximation. Following this, our general parametric solution is presented in section 4. Then in section 5, we prescribe a specific substrate material, and values for the vdW interaction energy between the substrate and graphene are obtained. In section 6 our results, supported by figures and numerical values, are presented and comparison is made with theoretical and experimental studies. In the final section, a brief summary is given and some future work regarding folding multi-layer graphene are discussed.

Model formulation
We now propose a general formulation to model the folding conformation of multi-layer graphene sheets supported on a substrate. The model is constructed based upon the mid-line of a multi-layer graphene stack. The mid-line represents a physical graphene layer when the total number of layers is odd (N = 2n + 1), or a notional layer when the total number of layers is even (N = 2n) for some integer n. As an example, figure 1 shows the geometry of folded 3-layer graphene stack located on a substrate. Further we consider a 2D problem due to the assumed translational symmetry in the fold direction. The solid curve in figure 1 represents the mid-line of the fold which consists of a curved region of a length L fold and a flat region of a length L flat = x 3 − x 2 .
The mid-line is also divided into three parts according to the sign of the line curvature and line gradient. The first part is denoted by C 0 which is from the point (x 0 , y 0 ) to the point (x 1 , y 1 ) with a strictly negative curvature. The second part, C 1 , is from the point (x 1 , y 1 ) to the point (x 2 , y 2 ) where the line curvature here is strictly positive. The flat region is C 2 from the point (x 2 , y 2 ) to the point (x 3 , y 3 ) with zero curvature. Also we use C to denote the concatenation of these three parts C 0 , C 1 , and C 2 .
The radius of curvature of the mid-line is denoted by with κ m denotes the associated line curvature given by where the curve is being parametrised in the curved region and dots here denote differentiation with respect to the natural parameter s. We assume that the layers maintain an equal distance 2δ g from one another and hence, the radius of curvature of the jth outer/inner curve is given by R j,o/i = R m ± νδ g , where j ä {1, K, n}, and ν = 2j when N is odd, or ν = 2j − 1 when N is even. Consequently, the corresponding jth outer/inner line curvature is given by Throughout this paper, δ g is used to denote the half spacing distance between any two closest graphene layers, and δ s is used to denote the spacing distance from the substrate to the first graphene layer.
The dominant energies here are the same as those considered when modelling the fold of a single-layer graphene sheet, namely the elastic bending energy and the vdW interaction energy [13]. Moreover, we model the elastic energy E e by integrating the sum of the squared curvature of each graphene layer over the total length of the curve C multiplied by its bending rigidity. In this study, we assume all graphene layers have the same bending rigidity. Thus, we may express the elastic energy as where γ is the bending rigidity of the multi-layer graphene stack, s is the arc length, and Generally, for each pair of inner and outer curves a distance νδ g away from the mid-line, we may combine their squared curvatures to simplify (1). Thus we may write Hence, (1) may be expressed in terms of κ m , the line curvature of the mid-line, as For either even or odd N, we deduce that , and 1 3 7 30 .

2
Therefore, we rewrite (3) as We model the vdW interaction energy between graphene and the substrate using a Heaviside unit step function u(x). In the folded region, we disregard the vdW interaction as an assumption of the modelling, but assume that it is the dominant contribution to the energy in the flat region. Thus the vdW interaction energy is modelled as where ò is a positive constant denotes the vdW interaction energy per unit area, and u(x) is a Heaviside unit step function. Noting that  = y 0 and ds=dx in C 2 . Thus, the total energy per unit length may be expressed as Furthermore, we apply an isoperimetric constraint on the total arc length of the curve C, that is Under this length constraint, a Lagrange multiplier λ is introduced into the total energy functional (5) which becomes In the above expression, the Lagrange multiplier λ plays no role in determining the mid-line that minimizes the functional E tot . So to simplify our model we let λ = ò in (6), and the functional to be minimized is now of the form subject to the following boundary conditions g , and the value of x 2 is to be determined from a natural boundary condition.

Variational calculus
Before considering the variation of (7), we first nondimensionalize it by introducing a scaling factor α such that x = αX and y = αY. In terms of the new variables, κ m and ds are now given byk k a = m m and dS = α ds, respectively. Making use of the new variables X and Y, the nondimensionalized energy functional is given by where S 2 is the nondimensional arc length at the point (X 2 , Y 2 ). Throughout this paper, we use F to denote the integrand part of (9), that is The variational problem for (9) leads to the pair of Euler-Lagrange equations where dots now denote differentiation with respect to the nondimensional arc length S.
which has no explicit dependence on X or Y. Thus, on integrating the pair of equations (11) with respect to S we obtain where β X and β Y are the integration constants. Since F has no explicit dependence on S, this provides another first integral given by which may be simplified by substituting β X and β Y from (12) as̈( where H is a constant. We may also derive the standard equation for the first variation ofĒ tot as Since X 2 , and therefore S 2 , is not prescribed at the end point S = S 2 , we require both β X = 0 and H = 0 due to the natural boundary condition. Thus (13) is now reduced tö̈( In the next three subsections we prescribe three different quantities for the integrand F based on different approximations ofk 2 , and for each approximation we apply (14) to derive the relative expression ofk m .

The one-term approximation
In this subsection we consider the simplest approximation of κ 2 with only the first term of (4) is included in the model, that is Under this approximation, the function F in (10) is now given by which after substituting the scaling factor α by g  N becomes Upon the substitutions ofF X ,F Y , and F from (15) into (14) we may derive whereupon we make the change of variables such that and primes denote derivatives with respect to θ. Therefore (16) may be rewritten as We now make the substitutionsˆ( sin , 1 7 m m which lead to an expression relating the nondimensional curvature given bŷ This expression is identical to that derived in [19], except that here we use a slightly different derivation and γ is now replaced by Nγ.

The two-term approximation
In this subsection we take into account an additional term in κ 2 with both the first and second terms of (4) are included in the model, that is We comment that for small δ g κ m = 1 this expression converges to the one-term approximation. Following the same steps as in section 3.1, we substitute a g =  N and rewrite the function F in (10) as Furthermore, we employ (14) to derive the following expression We now make the change of variables as in the previous subsection to obtain which may be rearranged to give ⎡ leads to the the nondimensional curvature of the mid-line, that is and by making the substitutions in (17), we deduce that

The three-term approximation
In this subsection we take into consideration a third additional term in κ 2 with all of the terms of (4) are included in the model, that is As before, we note that for small δ g κ m = 1 this expression approaches the two-term and one-term approximations. Again, the same substitution of α is used, and the function F in (10) is now given by where μ is given by (19), and After the substitution of the new integrand F into (14), we may compute the following expression Making the same change of variables as in the last two subsections we derive or equivalentlyˆ( which after some rearrangement reduces (22) to , and 2 5 5 1 5 1 . Now we apply Cardano's formula, and the solutions of (23) we seek are of the form w w w = + , 1 2 where 27 2 , and 4 27 2 . These solutions must satisfy (23) as well as the relation ω 1 ω 2 = − p/3. It may be shown numerically that the discriminant of the equation, D = − 27q 2 − 4p 3 , is always negative. Therefore we get three distinct roots, two are conjugate complex roots, and the third is a real root which may be written as Thus, using the substitutionˆ( ) we deduce that the line curvaturek m iŝ where p, μ, and ρ are parameters as defined earlier in this section, but using the substitutions in (17), q here is being a function of θ given by

Parametric solution
After determining the expression ofˆ( ) k q m for each approximation of κ 2 , we now use the two differential equations given in (17) to derive parametric solutions of the folded region of the curve C. Analytical solutions, for the one-term approximation, and numerical solutions, for the two-term and three-term approximations, are presented in the next two subsections.

The one-term approximation
The expression of line curvatureˆ( ) k q m given in (18), has been used by the present authors to derive the analytical solutions for the conformation of a rippled graphene sheet [19]. In the present paper, similar derivations to those in [19] are adopted to derive the analytical solutions for the folding conformation of supported multi-layer graphene stack. The substitution θ = 2f − π/2 is used to simplify our calculations, and the obtained solution is significantly driven by 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 modulus For C 0 , the part of the solution with negative curvature, the curve ranges between the starting point (X 0 , Y 0 ) at θ = θ 0 and the point of zero curvature (X 1 , Y 1 ) at θ = θ 1 . We note that θ 0 = − π which implies f 0 = − π/4, and the value of ( ) f =k sin 1 1 1 is determined by solvinĝ ( ) k f = 0 m 1 for f 1 . The solution for this part is given by  Y 1 ) wherek m changes sign. To satisfy this constraint, we repeat the same derivations as for the curve C 0 but use the opposite sign ofˆ( ) k f m . The solution for this part is given by 1 1 where f ä [f 1 , π/4].

The two-term and three-term approximations
In this subsection, we express our parametric solutions of the two-term and the three-term approximations in integral forms which are evaluated numerically due to the complexity of the corresponding expressions ofˆ( ) k q m given in (21) and (24). We comment that the line curvatureˆ( ) k q m changes sign at the point (X 1 , Y 1 ) which corresponds to θ = θ 1 . Thus by solvingˆ( ) k q = 0 m 1 for θ 1 , we may deduce that for each approximation. For C 0 , the part of the solution with negative curvature, we have θ ä [ − π, θ 1 ] and by taking into consideration the boundary conditions at the starting point θ 0 = − π, we deduce that For C 1 , the part of the solution with positive curvature, we seek a continuous curve at the point (X 1 , Y 1 ) wherek m changes sign. To satisfy this constraint, we repeat the same steps as for the curve C 0 but use the opposite sign of ( ) k q m . The solution for this part is given by sin d , Furthermore, we may obtain the solution of the jth outer/inner curve using where X(θ) and Y(θ) represent the parametric solutions of the mid-line, and ν is defined earlier in section 2. Multiplying these solutions by the scaling factor a g =  N recovers the dimensional solutions for the folded curve of C. Now the model behaviour is governed by material parameters related to the supporting substrate which are determined in the following section.

The supporting substrate
The silicon dioxide substrate (SiO 2 ) has been widely used as a supporting substrate in the fabrication of graphene-based electrical devices [21]. Also, the existing experimental data reported in the literature on folding multi-layer graphene sheets on SiO 2 [22], provides a good opportunity to validate our model. Based on these considerations, we specialize our solution to determine the folding behaviour of multi-layer graphene sheets supported on SiO 2 substrate. The Lennard-Jones (LJ) potential is used to account for the interaction energy between a pair of non-bonded atoms, which is written as where ξ and σ are known as the LJ potential parameters which are empirically determined, and j is the distance between atoms. The parameters A = 4ξσ 6 , and B = 4ξσ 12 are positive constants of the attraction and repulsion, respectively. For different number of graphene layers, we calculate the interaction energy ò between graphene and SiO 2 using the LJ potential which we approximate as where ò 2Ngs , ò Ngs , and ò Ng , denote the interaction energy between 2N graphene layers and the substrate, N graphene layers and the substrate, and N graphene layers, respectively. For the interaction energy between graphene layers, we apply a formula derived in Baowan et al [23], Ch. 3, that is where A gg and B gg are the graphene-graphene attractive and repulsive constants, respectively. Considering a three-dimensional SiO 2 substrate, we derive a formula for the interaction energy between graphene layers and the substrate by integrating (25) from 0 to ∞ with respect to z and replacing δ by δ + z. As the interactions between graphene layers and the SiO 2 substrate involve graphene to be interacting with both silicon and oxygen atoms, we deduce that with A gSi/gO and B gSi/gO being the graphene-silicon/graphene-oxygen attractive and repulsive constants, respectively. The LJ parameters,  [24]. In this study, η C = 0.3812 Å −2 is used to denote the surface density of carbon atoms on a sheet of graphene, where η Si = 0.02654 Å −3 and η O = 0.05308 Å −3 are used to denote the volume densities of silicon and oxygen atoms in the SiO 2 substrate, respectively.
Also in this calculation we take into account the change in the equilibrium distances as the number of layers changes. The values of δ g and δ s are determined such that (25) and (26) takes a minimum value, respectively. For example, the value of δ g = δ is chosen such that dò Ng /dδ = 0 is satisfied, as illustrated by figure 2. Table 1 summarizes the values we obtain for δ g , δ s , and ò for different number of graphene layers folded on a SiO 2 substrate. Substitution of the parameter values from table 1 into our solutions presented in section 4, leaves only one parameter still to be determined, namely β Y . Using the boundary condition y 2 = δ s + (3N − 1)δ g , we may calculate the value of β Y , and hence fully determine our solution.

Results
In this section, we first investigate the effects of different approximations of κ 2 on the folding conformation of multi-layer graphene located on a SiO 2 substrate. The values of bending rigidity γ = 1.0 eV, and number of layers N = 2 are adopted for this analysis. In figure 3, we present the predicted profile of the folded 2-layer graphene sheets on a SiO 2 substrate, for each considered approximation of κ 2 . The analysis reveals remarkable  differences between these approximations in the hump height of the fold, that is the measured height from the point where the folding edge becomes flat to the maximum height obtained. Using the one-term approximation, we obtain a hump height of h = 0.93 Å which is significantly smaller than the hump heights of h = 5.23 Å, and h = 5.68 Å which are obtained using the two-term and the three-term approximations, respectively. We note that the one-term approximation in this work is not shown to be consistent with previous experimental measurements [22], and thus is not used in further comparison. However, we comment that the one-term approximation is largely equivalent to the model used by Meng et al [25] where the folding conformation of unsupported multi-layer graphene is modelled. The folding profile of 8-layer graphene is accurately predicted in [25], using both the finite-deformation model and molecular dynamics simulations. The difference in predicting the folding profile may be justified by the value of the vdW interaction energy used in each model. Since this model is applied to supported multi-layer graphene, the vdW interactions between graphene layers and the substrate are relevant here. Further, we compare our solutions to experimental measurements of hump heights. The calculated hump heights of 2-to 6-layer graphene sheets folded on a SiO 2 substrate are shown to be consistent with previous results [22], and a representative plot for the comparison is shown in figure 4. As the number of layers N increases, our model predicts the hump height more accurately as more terms are included in the curvature squared approximation κ 2 . However, for N 6, it appears that even the three-term approximation would be insufficient to accurately predict hump heights. For reasons of space, selected range of folding profiles are presented in figure 5 for different number of layers N. Although we do not take into consideration any defects in the graphene structure, our model shares a number of similarities with the experimental data reported by Chen et al [22].

Summary
In this paper, we develop a mathematical model to determine the folding conformation of multi-layer graphene sheets located on a substrate. Due to the assumed translational symmetry of the fold, a 2D problem is proposed. The model is nondimensionalized, and the calculus of variations is used to minimize the energy functional and determine the shape of the fold. Solutions are derived for different approximations of κ 2 based upon the midline, which is introduced as a notional curve when the number of layers is even. We consider SiO 2 as the supporting substrate, and the vdW interaction energy between graphene layers and the substrate as well as the equilibrium spacing distances are obtained using LJ potential and presented in table 1. The boundary condition on y at the end point s 2 is used to determine the value of the unknown constant β Y that arises when integrating the Euler-Lagrange equations.
The folding conformation of multi-layer graphene located on a SiO 2 substrate is investigated using three different approximations of κ 2 . The one-term approximation yields much smaller hump height than the other two approximations. Also, hump heights obtained from our model are adopted to compare with the experimental data reported by Chen et al [22] for different number of layers. Representative plots for predicted folding profiles and hump heights comparison are shown in figures 3, 4, 5. Whilst we apply this model to supported multi-layer graphene and account for the interaction energy between graphene layers and the substrate, our model is shown to be consistent with experimental data for a maximum of 6 layers. We mention that earlier models have looked at the folding of unsupported multi-layer graphene [15,25], where the folding profile of up to 8-layer graphene is predicted using theoretical models and simulations. As we considered up to three terms in the expression of κ 2 , a possible future work may be adding extra terms in order to investigate the folding conformation of supported multi-layer graphene with more than 6 layers. However, simply extending the Taylor series would certainly complicate the derivation of the corresponding expression for the line