Variational model for a rippled graphene sheet

The calculus of variations is utilised to study the behaviour of a rippled graphene sheet supported on a metal substrate. We propose a model that is underpinned by two key parameters, the bending rigidity of graphene γ, and the van der Waals interaction strength ξ. Three cases are considered, each of which addresses a specific configuration of a rippled graphene sheet located on a flat substrate. The transitional case assumes that both the graphene sheet length and substrate length are constrained. The substrate constrained case assumes only the substrate has a constrained length. Finally, the graphene constrained case assumes only the length of the graphene sheet is constrained. Numerical results are presented for each case, and the interpretation of these results demonstrates a continuous relationship between the total energy per unit length and the substrate length, that incorporates all three configurations. The present model is in excellent agreement with earlier results of molecular dynamics (MD) simulations in predicting the profiles of graphene ripples.


Introduction
Graphene is a two-dimensional sheet of carbon atoms bonded to each other in a planar hexagonal array. This two-dimensional structure endows graphene with many useful electronic and mechanical properties. [1][2][3][4][5][6][7] These properties mean graphene is a highly promising material for constructing nanoelectromechanical systems. 8,9 Graphene is hypothesised to be a biocompatible material since its bending stiffness is comparable with that of the lipid bilayers of the biological cells, 10 and it has many other applications in gas separation, 11 biomedicine, 12 nitrogen reduction reaction, 13 and metal-ion batteries. 14 For example, popgraphene is a benecial anode material in lithium-ion batteries with fast charge and discharge rates. 15 However, graphene does not always remain planar. For example, ripples are observed in suspended graphene sheets, 16 as well as graphene on a substrate. 17 Experimental studies nd that the electronic properties of graphene can be affected by the range and height of these ripples. 18,19 Gui et al. 20 use rst-principles calculations to predict the electronic properties of a rippled graphene sheet where a band gap opening is observed in the rippled graphene. They evaluate a direct band gap at 0.93 eV which indicates rippled graphene may be a highly tunable semiconductor. They also report that the band gap increases proportionally with ripple amplitude, which is the maximum distance between the graphene sheet and the substrate.
The pattern of ripples in graphene may be affected by a variety of factors, including temperature, substrate material, and the size of the graphene sheet. For example, suspended graphene remains at with no ripples when the temperature is close to absolute zero, but ripples appear as the temperature increases. 21 When a graphene sheet is placed on a substrate the amplitude of the ripples depends on the substrate properties such as roughness and interfacial van der Waals interactions. 17 Moreover, the overall size of the graphene sheet has a signicant impact on the amplitude of ripples, with the ripple amplitude increasing proportionally with the size of the sheet. 22 The bending rigidity of graphene plays an important role in the conformation of ripples. Using density functional theory calculations, Wei et al. 10 evaluate the bending rigidity of a single-layer graphene at 1.44 eV. They compare this result to earlier obtained values for the bending rigidity of graphene which range from 0.80 to 1.60 eV. To cover this range of values, we will consider various bending rigidities for the rippled graphene in the present paper.
The energies we consider when modelling ripples in graphene are the same as those considered when modelling a fold of graphene. Cox et al. 23 develop a model which describes a graphene sheet folded onto itself where two energies are taken into account. The rst is the elastic energy that arises when bending a graphene sheet which is represented mathematically by integrating the square of the line curvature over the length of the total curve and scaled by the bending rigidity of graphene. The second energy is the van der Waals (vdW) interaction between the at section of the graphene layers, prescribed as the integral of a Heaviside unit step function multiplied by the vdW energy per unit area. An extension of the fold model was further developed by Dyer et al. 24 where either a carbon nanotube of circular or elliptical cross section is encapsulated in the graphene fold. A comparison between MD simulations and the mathematical model shows excellent agreement between the two approaches in determining the fold shape.
In the following section, we formulate a general mathematical model of a single graphene ripple on a substrate. In Section 3 we non-dimensionalise this model and employ the calculus of variations assuming a xed length (isoperimetric) constraint of both the rippled graphene sheet and the substrate. In Section 4 we relax the isoperimetric constraint on the sheet length, but maintain a xed substrate length. The substrate length is varied in Section 5 and an isoperimetric constraint is applied to only the sheet length. The solutions and associated numerical details are provided with each case. The main results and the relationships between the three cases considered here are described in Section 6. Discussion and some concluding remarks are made in Section 7.

Model formulation
We propose a general formulation to determine the conformation of a rippled sheet of graphene on a substrate for each of the three cases mentioned above. We consider two different metal substrates: Cu(111) and Ni(111). The geometry of the rippled graphene sheet is shown in Fig. 1. Here we assume a translational symmetry in the z-direction which reduces the problem to two dimensions. Also, the reective symmetry about the yaxis allows us to consider only solutions in the right half plane. We further divide the solution curve into three sections. The rst section C 1 is the curve from the point (0, h) to the point (x 0 , y 0 ) where the line curvature is strictly negative. The second section C 2 is from (x 0 , y 0 ) to (x 1 , d) where the line curvature is strictly positive. The third section C 3 is the horizontal line from (x 1 , d) to (x 2 , d) where the line curvature is zero. The concatenation of these sections is denoted by C ¼ C 1 + C 2 + C 3 . In these position vectors we use h to denote the ripple height, and d to denote the separation distance between the substrate and the at section of the graphene sheet.
With reference to the model developed by Cox et al., 23 the dominant energies for our model are the elastic bending energy and the vdW interaction energy. The elastic energy E e is modelled with where g is the bending rigidity of graphene, s is the arc length, and k is the line curvature of y ¼ y(x) which is given by As mentioned in Section 1, the bending rigidity of graphene reported by various authors ranges from 0.8 to 1.6 eV as presented in ref. 10, and for the purpose of comparison we consider a linear sample from this range for g, namely, 0.8, 1.0, 1.2, 1.4, and 1.6 eV. We model the vdW interaction energy between the graphene and the substrate as where x denotes the vdW energy of graphene-substrate interactions per unit area, and u(x À x 1 ) is a Heaviside unit step function. Our model assumes that the elastic energy dominates in sections C 1 and C 2 , and the vdW interaction dominates in C 3 . We comment that the curvature and therefore the elastic energy vanishes in C 3 , but discarding the vdW interaction in sections C 1 and C 2 is an approximation we make in the modelling. Therefore, the total energy per unit length may be expressed by subject to the boundary conditions where at the endpoint x ¼ x 1 the value of x 1 is not prescribed and we have a natural boundary condition on x. Here dots denote differentiation with respect to x.
For the purpose of comparison we prescribe a particular value of h for each case of a rippled graphene sheet supported on a different substrate, and the reason for adopting these particular values will be explained in Section 4. The prescribed ripple height, h, and empirical values for the separation distance between graphene and the substrate, d, and the vdW energy of graphene-substrate interactions per unit area, x, corresponding to each substrate are shown in Table 1. Now we apply this model to consider the three cases for the length of graphene sheet placed on a substrate, namely, the transitional case, the substrate constrained case, and the graphene constrained case. Fig. 2 illustrates the geometry of each case.

The transitional case
We now apply the calculus of variations to determine the conformation of a rippled graphene sheet and derive a solution   (1) is minimised. We impose an isoperimetric constraint on the total arc length of C, given by ð C ds ¼ L: Moreover we assume that x 2 is xed, but x 1 varies (to be determined from a natural boundary condition). Therefore the new functional we wish to minimise is of the form subject to the boundary conditions (2), where l is a Lagrange multiplier corresponding to the isoperimetric constraint. We nondimensionalise (3) by dening X ¼ x/a and Y ¼ y/a where a is a scaling factor, and thereforeÿ ¼ Y 00 /a and _ y ¼ Y 0 , where primes denote derivatives with respect to the nondimensional X coordinate. The substitution of these new variables into eqn (3) yields þ aðl À xÞðX 2 À X 1 Þ: Furthermore we subtract the xed energy at X 2 , and divide both sides of eqn (4) by ax. Thus we derive the nondimensionalised energy functional aŝ We let a ¼ ffiffiffiffiffiffiffi ffi g=x p and m ¼ l/x so that which is the functional we wish to minimise subject to the boundary condition Y(0) ¼ h/a and a natural boundary condition applies at X ¼ X 1 . To simplify the following calculations, we use f to denote the integrand part of (6), that is Extremals of (6) are given by Euler-Lagrange equation and since (7) does not depend on Y explicitly, then on integrating (8) with respect to X, we obtain the rst integral where b is a constant. Further, the integrand has no explicit dependence on X, this provides the additional rst integral where H is a constant. On considering the corresponding second-order variational problem, we may derive the standard equation for the rst variation ofÊ tot as where Here the natural boundary condition, which applies when the Xcoordinate at the end point X ¼ X 1 is not prescribed, requires H ¼ (1 À m). Therefore by combining eqn (9) and (10), we obtain Which aer the substitution of (7) leads tô which is an equation relating the nondimensional curvaturek and the rst derivative of Y. In order to obtain the parametric solutions, we make the substitution of Y 0 ¼ tan q which yieldŝ The above expression may be written as followŝ whereupon we dene a new parameter j, such that This leads to the succinct formulâ Now, using the fact thatk ¼ cos(q)dq/dX ¼ sin(q)dq/dY, we derive two rst order differential equations, namely To facilitate further integration we change the parametric variable of our equations from q to f such that q ¼ 2f À j, thus eqn (12a) and (12b) take the following forms 2ðcos j cos 2f þ sin j sin 2fÞ ðm þ ð1 À mÞsec j À 2ð1 À mÞsec j sin 2 fÞ 1=2 # ; (13a) 2ðcos j sin 2f À sin j cos 2fÞ ðm þ ð1 À mÞsec j À 2ð1 À mÞsec j sin 2 fÞ 1=2 # ; (13b) and we now obtain the general solutions by integrating eqn (13a) and (13b). This yields with the parameters A ¼ 2ðm þ ð1 À mÞ sec jÞ 1=2 ð1 À mÞ sec j ; B ¼ m m þ ð1 À mÞ sec j ; ; where F(f,k) and E(f,k) represent the incomplete elliptic integrals of the rst and second kind, respectively, with elliptic modulus k, and c 1 , c 2 are arbitrary constants of integration. We now dene two functions which play a signicant role in determining our parametric solutions, namely The line curvature changes sign at f 0 ¼ sin À1 (1/k) corresponds to (X 0 , Y 0 ) which is the boundary point between the two curves C 1 and C 2 , which has coordinates Moreover, the solutions are shown to have a rotational symmetry about this critical point (X 0 , Y 0 ) where f varies over the range [j/2, f 0 ]. Thus our solutions may be written in terms of X 0 and Y 0 as Now, using the boundary condition Y C 2 (j/2) ¼ d/a, we may obtain d À h 2a ¼ À½sin j g 1 ðj=2Þ þ cos j g 2 ðj=2Þ: Hence we may rewrite eqn (15b) as which introduces the Y component of the critical point as the midpoint of the ripple amplitude and the at section of the graphene sheet. We also comment that x ¼ aX and y ¼ aY, and so multiplying these solutions by the scaling factor a ¼ ffiffiffiffiffiffiffi ffi g=x p recovers the re-dimensionalised solutions.

Numerical results
We are le with two parameters to determine, namely j and m.
We may nd their values by solving the following system of equations numerically y C 2 ðj=2Þ ¼ d; and where d takes a specic empirically determined value for each substrate as given in Table 1, and L is the assumed xed arc length of C. Aer nding the unknowns j and m, the solution is fully determined and representative plots are shown in Fig. 3. At this point, we may also give the total half arc length of the ripple L rip by integrating ds over the curve This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 16016-16026 | 16019 Paper RSC Advances Again using the substitution of q ¼ 2f À j we deduce that where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m þ ð1 À mÞsec j p . Hence, the total half arc length of the ripple is given by We may also calculate the elastic energy by integrating the square of the line curvature over the total length of C 1 + C 2 , as follows and therefore, the elastic energy is given by E e ¼ 4r ffiffiffiffiffi gx p ½Eðf 0 ; kÞ À Eðj=2; kÞ: Hence, the total energy per unit length for this model is given by ffiffiffiffiffi gx p ½Eðf 0 ; kÞ À Eðj=2; kÞ À xðx 2 À x 1 Þ: Numerical values for the x-component of the critical point x 0 and the total half arc length of the ripple L rip are presented in Table 2, for a range of values of the bending rigidity g. The vdW interaction x takes a specic empirically derived value for each substrate as given in Table 1.

The substrate constrained case
We now consider the special case, shown in Fig. 2(b), by removing the assumed isoperimetric constraint on the total arc length of C to allow graphene to overhang the substrate. Therefore we discard the variation of x 2 and substitute l ¼ 0, and a ¼ ffiffiffiffiffiffiffi ffi g=x p into eqn (5) to obtain the functional that we wish to minimise, that iŝ We note that this expression is identical to (6) with m ¼ 0, and so by the same reasoning as in Section 3 we obtain k ¼ AE[sec j cos(j + q)] 1/2 .
Consequently, we redene some of our parameters so that . Again the line  curvaturek changes sign at f 0 ¼ sin À1 (1/k) ¼ p/4 which corresponds to the point (X 0 , Y 0 ) and has similar coordinates to that in (15a) and (15b). Taking the redened parameters into account, these coordinates are Similarly, with the rotational symmetry of the solutions about the point (X 0 , Y 0 ), and the variation of f over [j/2, p/4], the curves C 1 and C 2 still can be written as Using the same derivation as in Section 3, the Y component of the critical point in this case may be rewritten as As before, multiplying the solutions above by the scaling factor a ¼ ffiffiffiffiffiffiffi ffi g=x p recovers the dimensional solutions for the curves C 1 and C 2 .

Numerical results
In this case we have only one parameter to determine, namely j, and using the boundary condition Y C 2 (j/2) ¼ d, the value of j may be determined numerically. Aer determining j, the solution is fully determined and representative plots are shown in Fig. 4. We comment here that the gradient of the graph for g ¼ 0.8 in the neighbourhood of the critical point becomes innite when the ripple amplitude h z 8.0 d for the Cu(111) substrate, and h z 6.7 d for the Ni(111) substrate where d takes a specic empirical value for each substrate as presented in Table 1. We note that by increasing the ripple height the graphene sheet starts to bend over itself which is the transition point from the rippled state to the wrinkled state. Also, higher bending rigidities require higher ripple amplitudes for the maximum gradient to approach innity. However, we x these values of h for each corresponding substrate in all plots to enable comparison between various values of the bending rigidity g.
The total half arc length of the ripple L rip and the total energy E tot for this case can be calculated from eqn (17) and (18), respectively, with the new parameters m ¼ 0 and f 0 ¼ p/4. Thus, the total half arc length of the ripple is given by and the total energy per unit length is given by ; kÞ À Eðj=2; kÞ À xðx 2 À x 1 Þ: Numerical values for the x-component of the critical point x 0 and the total half arc length of the ripple are shown in Table 3,  (111) substrate, for various bending rigidity g, in the substrate constrained case. Table 3 The x-component of the critical point x 0 and the ripple half arc length L rip for various bending rigidities g, in the substrate constrained case g (eV) Cu (111) substrate Ni(111) substrate for a range of values for the bending rigidity g. Again the vdW interaction strength x takes a specic empirically derived value for each substrate as given in Table 1.

The graphene constrained case
In this case, as shown in Fig. 2(c), in addition to the assumed isoperimetric constraint on the total arc length of C, we take into account the variations of x 2 , and assume it varies proportionally with x 1 , that is Following the nondimensionalisation described in Section 3, the total energy per unit length for this case iŝ where ds ¼ a(1 + Y 02 ) 1/2 dX. We now divide both sides by xa and set a ¼ ffiffiffiffiffiffiffi ffi g=x p . Further, the aim is to determine the shape of Y ¼ Y(X) which minimisesÊ tot , and here l has no impact on the solution, so we choose l ¼ x to simplify the calculation. The nal form of the functional to be minimised iŝ subject to the boundary conditions provided in Section 2.
Comparing (21) and (6), we nd that m ¼ 1. In a similar way as in the previous two cases, the use of calculus of variations then provides thatk Similarly, by comparing (22) with the equivalent expression for the transitional case (11), we note that in (22) j ¼ p/2 and b acts in place of the term ð1 À mÞsec j ¼ as m approaches 1. Therefore, we redene the parameters A, B, and k in (14a) and (14b) as : The critical point (X 0 , Y 0 ) represents the point of the rotational symmetry which again corresponds to f 0 ¼ sin À1 (1/k), and the coordinates of this point are given by Similarly, the solutions on the curves C 1 and C 2 where f[ f 0 , p/4] are Now, we use the boundary condition Y C 2 (p/4) ¼ d/a to rewrite the Y component of the critical point as  (111) substrate, for various bending rigidity g, in the graphene constrained case.
The dimensional solutions for the curves C 1 and C 2 may be obtained by scaling these solution by a ¼ ffiffiffiffiffiffiffi ffi g=x p .

Numerical results
Here again we are le with only one parameter to determine, namely b, and using the boundary condition Y C 2 (p/4) ¼ d, the value of b may be determined numerically. Aer nding b, the solution is fully determined and representative plots are shown in Fig. 5. The total half arc length of the ripple L rip for this case can be evaluated from (17) with j ¼ p/2 and r ¼ ffiffiffiffiffiffiffiffiffiffiffi Similarly, the total energy per unit length E tot for this case is obtained from (18), which yields E tot ¼ 4r ffiffiffiffiffi gx p ½Eðf 0 ; kÞ À Eðp=4; kÞ À xðx 2 À x 1 Þ: Numerical values for the x-component of the critical point x 0 and the total half arc length of the ripple L rip are presented in Table 4 for a range of values of the bending rigidity g. As before, the vdW interaction x takes a specic empirically derived value for each substrate as given in Table 1.

Results
For the purpose of comparison, we assume a xed half length L ¼ 700Å for the rippled graphene sheet on a at substrate. The goal is to obtain a relationship between the variation of the substrate length and the total energy per unit length for the rippled graphene. The substrate half length L sub is divided into three regimes that correspond to the point x 2 , where the graphene edge ends, which can be calculated by Here we comment that x 2s < x 2g where x 2s and x 2g correspond to the substrate constrained case and the graphene constrained case, respectively.
The substrate constrained case was formulated to describe the regime where L sub < x 2s , provided that as L sub gets closer to x 2s the total energy per unit length E tot decreases due to the vdW interactions that occur when L sub increases. The graphene constrained case can be adapted to address the regime when L sub > x 2g , and here the total energy per unit length E tot remains constant since there are no more interactions that are involved for this section. Finally, the transitional case is used to evaluate the total energy per unit length E tot when L sub z x 2 . Multiple values for x 2t were chosen from the interval [x 2s , x 2g ] to ll the gap between the substrate constrained case and the graphene constrained case where x 2t denotes the location of the graphene edge in the transitional case. Fig. 6 illustrates the relationship between varying the substrate length and total energy per unit length included the three cases discussed above.
A graphene sheet of half length L ¼ 700Å will remain at when it is placed on a at substrate of the same half length. The substrate constrained case can be adapted to investigate the impact of reducing the substrate length on the graphene sheet. Fig. 7 shows that ripples form in this scenario, and that there is an increase in the ripple amplitude as the substrate length decreases. Graphene of lower bending rigidity reaches the Table 4 The x-component of the critical point x 0 and the ripple half arc length L rip for various bending rigidities g, in the graphene constrained case g (eV) Cu (111)     translational symmetry along the ripple reduces our problem to two dimensions, and the reective symmetry about the y-axis allows us to consider only solutions in the upper right quadrant. We employ variational calculus to minimise the elastic energy arising from the curvature squared and maximising the vdW interaction energy between the at section of the graphene and the metal substrate. We account for the length of the substrate by considering three different cases. The rst addresses the case when the edges of the rippled graphene sheet and the substrate coincide. The second case is when the edge of the graphene sheet overhangs the substrate edge, and the last case applies to a rippled graphene sheet for which the edge of the sheet does not extend to the substrate edge. We assume a xed half length of the graphene sheet, and consider all three cases to obtain a continuous relationship between the total energy per unit length of the graphene and the length of the substrate as shown in Fig. 6. The substrate constrained case is used to illustrate the effect of reducing the substrate length on the ripple height and the total energy of the graphene as shown in Fig. 7 and 8. Using the transitional case, our model is shown to agree with the results of earlier MD simulations of graphene ripples (or wrinkles) on a Cu(111) substrate.
Future research may address the phenomenon of ripple formation when the substrate shrinks by taking into account additional physical effects. For instance the model maybe extended to account for the change in the vdW interaction strength as the substrate length and density change. This will introduce x as a function of x and lead to an explicit dependence on x in the integrand part of (6), which would undoubtedly complicate the corresponding Euler-Lagrange equation.

Conflicts of interest
There are no conicts to declare.