A variational model for conformation of graphene wrinkles formed on a shrinking solid metal substrate

Chemical vapor deposition is a popular technique for producing high-quality graphene sheets on a substrate. However, the cooling process causes the graphene sheet to experience a strain-induced, out-of-plane buckling. These wrinkles structures can have undesirable effects on the properties of the graphene sheet. We construct a pair of models to analyse the conformation structure of these wrinkles. An arch-shaped wrinkle is first modelled then expanded to incorporate self-adhesion between the wrinkle edges. Variational techniques are employed on both models to determine the optimal conformation for graphene supported on Cu and Ni substrates. We find these models predict a similar structure to experimental analysis of graphene wrinkles on these solid metal substrates.


Introduction
Graphene is often visualised as a perfectly flat material. Contrary to this assumption, graphene sheets typically exhibit varying surface corrugations. Structural defects result in bonding asymmetry that is relaxed by an out-ofplane deformation. Even relatively pristine free-standing graphene is still affected by thermal vibrations resulting in ripple formation and self-folding at the edges. In particular, ripples strongly affect many properties of graphene, including electronic properties, chemical properties and electrical performance [1][2][3][4]. Park et al [4] find that with a high cooling rate in chemical vapor deposition (CVD), graphene sheets have reduced ripple density and heights, which lead to enhanced electrical characteristics, such as higher electron/hole transports and reduced sheet resistance. Commonly, graphene grown on a substrate is subject to wrinkle formation from in-plane stresses. Wrinkled structures can also be formed by the overlapping of graphene islands during production. Additionally, high-densities of wrinkles are formed by compressive strain induced by thermal effects during the cooling process [5][6][7]. Not only during the production process of graphene, wrinkles are also observed after graphene is transferred onto other substrate for applications [8].
Deformations of graphene greatly impact the properties of the sheet. Amongst other effects, graphene wrinkles cause anisotropic electron mobility, potential band-gap opening, locally-modified chemical reactivity and a reduction in both mechanical strength and thermal conductivity. The impact of these effects largely depends on the desired applications. Tuneable wrinkles enable functionalisation of the sheet and have potential application in energy storage and selective bandgap engineering [9]. However, wrinkles are generally regarded as structural defects with usually deleterious effects.
We examine wrinkled structures induced during CVD. Wrinkle formation is inherent to this process. Subsequently, understanding the nature of these wrinkles remains an important research objective [10]. Chemical vapor deposition is a popularly used technique for producing high-quality, large-area graphene sheets. This method involves the bottom-up process of growing of graphene on a crystallographic substrate. Methane (or other carbon-rich feedstock) gas is introduced into a reaction chamber containing a substrate. Graphene produced on transition metal substrates can be categorised as weak binding from physisorption on Cu, Al, Pt, Ag and Au as well as strong binding from chemisorption occurring with Ni, Co, Ru, Pd and Ti interaction. We note that liquid metals and alloys, such as Ga, In, Sn, In-Cu, Sn-Ni and Sn-Ag-Cu and molten Cu and Ni are also used to grow graphene in CVD [11][12][13]. As mentioned in [11], the surface fluctuation and graphene wrinkles are also observed on liquid metals which are caused by different thermal expansion coefficients of liquid metals and graphene. Controlling the difference of these coefficients may lead to the ability to tune the morphology and graphene wrinkles and ripples on liquid metal surfaces [7,11].
In this paper, we are primarily interested in solid metal substrates. Although the model applies to general solid metal substrates, we perform our calculations for Cu(111) and Ni(111) substrates to compare our mathematical prediction with existing experimental and computational studies. In CVD, high temperatures within the chamber cause the feedstock to decompose and the resultant carbon atoms attach to the substrate. The solubility of carbon into copper is insignificant, resulting in a self-terminating process producing uniform mono-layer graphene over most of the surface [14,15]. As the graphene layer and the substrate cool, the lattice of both materials contract. The substrate shrinks notably, while the sheet expands slightly with the end result of compressing the graphene mono-layer laterally. In addition, the van der Waals (vdW) adhesion between the graphene sheet and substrate is weak and highly prone to interlayer slipping and delamination. Therefore, sufficient compression strain is relieved through out-of-plane buckling and localised wrinkle formation occurs in the graphene layer as seen in figure 1 [16,17]. Wrinkles tend to form around defects and defect lines in the substrate [17][18][19]. These wrinkles and other folded structures have been identified through scanning electron microscopy (SEM) and atomic force microscopy (AFM). Furthermore, graphene growth via CVD on copper substrate using methane is strongly affected by hydrogen (H 2 ) [20]. In particular, the pressure and flow of H 2 during CVD can control ripple density on graphene sheet and trapped H 2 molecules in the high-curvature convex regions of graphene ripple can influence the shape of wrinkles [21][22][23]. However, this paper does not consider the effect trapped H 2 molecules in formulating our mathematical model for wrinkles. Nevertheless, the incorporation H 2 effect in formulating mathematical model for wrinkle structures certainly opens an avenue for future research in this area.

Modelling approach
The calculus of variations has previously been employed for formulating continuous models for the joining of carbon nanostructures, which are taken to deform in perfect elasticity [25]. This approach has also proven successful at determining optimal graphene folding conformations, giving results consistent with molecular simulations and experimental results [26,27]. Recently, this approach is also used to determine conformations of rippled graphene sheet based on the constraints of the length of graphene sheet and the length of metal substrate [28]. In this paper, similar techniques to those of [28] are adopted to predict the structure of graphene wrinkles formed on substrate. Specifically, we examine the buckling resulting from lateral strain applied to the sheet. Wrinkle configurations are taken to be induced from the energy balance between the bending rigidity of graphene and the van der Waals adhesion between the sheet and substrate. We comment that this approach is to average across all possible conformations, and so the model solution can be considered a likely indicative solution. Individual cases may vary due to stochastic effects.
In the following subsections, we consider mathematical models for two graphene wrinkle profiles as shown in figure 2. We first examine a small delamination (Model I), characterised by an 'arch' shape. Under increasing strain, the sides of the wrinkle are found to approach one other. Subsequently, we examine a second configuration (Model II) where the sides of the wrinkle self-adhere through van der Waals interactions. Such a structure has been previously speculated in the literature [2].

Model I-Arch wrinkle structure
We construct a model to determine the optimal structure of a graphene wrinkle induced on a substrate through compressive force. We assume that a pristine sheet of graphene is initially attached to a perfectly flat substrate through van der Waals adhesion forces. This substrate contracts under negative thermal expansion inducing work upon the sheet. Under sufficient work, the graphene sheet experiences an out-of-plane buckling forming a wrinkle structure. Transmission electron microscope (TEM) imaging suggests a one-dimensional nature to these wrinkles [18]. Subsequently we take the graphene wrinkle to be invariant in the z-direction, reducing our problem to a two-dimensional curve as shown in figure 2(a). Furthermore, we assume reflective symmetry in the y-axis and therefore solely consider the solution curve the first quadrant.
The dominant energies are taken to be comprised of the elastic bending energy of the graphene sheet and the van der Waals interaction energy between the sheet and the substrate. The sheet is initially situated at the equilibrium distance δ 1 from the substrate. We additionally introduce a fixed lateral compression 2ℓ to the sheet.
The graphene is modelled in two parts. Firstly the curve  0 represents the free sections of the wrinkle, while the line from (x 1 , δ 1 ) to (x 2 , δ 1 ) represents the horizontal section of the graphene which remains adhered to the substrate. The curves denoting the wrinkle  0 is bounded by the points (0, h) and (x 1 , δ 1 ). The height h and width x 1 of the wrinkle are not prescribed a priori.
The elastic bending energy per unit length, E e is proportional to the integral of the squared curvature over the curve  0 and can therefore be modelled as where κ is the line curvature and γ is a constant denoting the bending rigidity of graphene. We model the van der Waals energy per unit length, E v , as = - where ò 1 denotes the van der Waals energy per area length. The total potential energy per unit length, E, after undergoing a gauge transformation of  x 1 2 and rescaling by γ/α 1 , yields the nondimensionalised functional Ê given by ò a g a g a k a = 1 , which is the characteristic length for the problem. We must incorporate a constraint related to the lateral compression forming the wrinkle. First, the graphene sheet is taken to be inextensible forcing the total arc length of the curve to be fixed at a length L. That is, Determining the minimal energy conformation is equivalent to minimising the functional J(y). We denote the integrand in This problem is of the form of a second order calculus of variations problem with the corresponding Euler-Lagrange equation, Since the endpoints (0, h) and (x 1 , δ 1 ) may vary in the y-and x-directions, respectively, we prescribe = Since F is independent of both x and y, both H and p are conserved quantities. Thus we deduce from the natural boundary conditions that p ≡ 0 and for extremals l a = -H 1 1 ( ) , over the whole curve  1 . The Euler-Lagrange equations therefore simplify to Upon substitution and some routine calculation, (5) may be rearranged to find, where the conformation of the curve is now expressed by its curvature. We therefore make the substitution, q ¢ = y tan , converting this problem into the standard parametric equations for a plane curve. Here, θ=θ(s) represents the tangential angle to the curve where s is the position along the arc length of the curve. We note that c o s , s i n .
Applying the chain rule, the curvature (6) can be expressed as With some rearrangement, we construct the parametric equations cos  1  cos  ,  sin  1 cos .
For ease of solving, we make the substitution θ=2f which leads to where we have simplified using the identity where d is a constant of integration. We now find a solution for x(f). Reorganising (7), we have and integrating yields where c is a constant of integration and F(f, k) and E(f, k) are the incomplete elliptic integrals of the first and second kind, respectively. The constants are determined from the boundary conditions. We first separate the parametric equations based on the sign of the curvature. The critical point of zero curvature corresponds to f=±f 0 , where f =k sin 1 0 1 ( ), and then endpoints where ¢ = y 0 corresponds to f=0. To construct the full extremal curve we begin by choosing the negative curvature solution and define the functions where the parameter f f Î -, 0 . The positive curvature solution has a rotational symmetry. To satisfy the endpoint conditions of the curve  0 starting at (x 1 , δ 1 ) and finishing at (0, h) we define parametric solutions for the two halves as where the subscript - 0 relates to the upper part of the curve where curvature is negative and the subscript +  0 relates to the lower part of the curve where curvature is positive. From these functions we may determine the wrinkle half-width x 1 and height h as given by Further recalling that the curvature κ=dθ/ds, and that q k f k = = ds d d 2 and substitution into the constraint equation (2) yields We also note that the elliptic integrals with argument f 0 can be easily expressed as Gaussian hypergeometric Using the expression for curvature and details of the length constraint we may write which again may be written with hypergeomtric functions as 1 2, 1 2; 1; 1 1 2, 1 2; 2; 1 , 12 Finally, we note that for larger choice of compression value, the wrinkle conformation forms a s-shaped or 'mushroom' structure [29,30]. For even larger compression values, a new model may be required, which we develop in the next subsection.

Model II-Self-adhered wrinkle structure
The previous subsection (Model I) captures the conformation of a small wrinkle induced upon an adhered graphene sheet. We note that a sufficiently narrow wrinkle could self-adhere to produce a more stable energyoptimal structure [2]. As such, we expand the previous formulation to account for this potential structural evolution. The curve representing the cross-section of this self-adhered wrinkle is shown in figure 2(b). The total profile comprises four parts: two curves labeled as  1 and  2 where elastic forces dominate; and two lines (one vertical at x=δ 0 for y 1 <y<y 2 , and one horizontal at y=δ 1 for x 1 <x<x 2 ) where van der Waals forces dominate. For this model we have the following total energy per unit length given by ò ò where γ is the flexural rigidity of graphene, ò 0 is the van der Waals energy per unit area for graphene-graphene interactions, and ò 1 is the van der Waals energy per unit area for graphene-substrate interactions. We now decompose this function into components as follows ò ò and we can then specify two solution curves, one for negative curvature --  x y , where for the negative curvature solution p f y -  . We now determine the energy contribution from the curve  1 . Recall that the elastic energy is given by The length of  1 is given by

Model II-Extremal for E 2
We recall that E 2 is given by and the solution is given by

Geometric parameters and total energy for model II
Apart from the material constants of δ 0 , δ 1 , γ, ò 0 and ò 1 we assume that the two geometric parameters that are specified a priori are the graphene half-length L and lateral compression half-length ℓ. From these two parameters we may immediately specify the graphene edge location, which is given by We also note that as already stated in the previous subsection we have determined the width and height of the base of the wrinkle which are given by x y y 2 , and 4 .
We also note that for Model II to be physical we require < y y 1 2 and x 1 <x 2 which may be expressed with the following pair of inequalities The total energy per unit length of the the conformation for Model II is given by

Results
Here, we construct a pair of models for graphene wrinkle conformation as detailed in section2 for Models I and II. Based on [31], we use the value of the bending rigidity, γ in the range 0.8-1.6eV. We study the most commonly employed substrates, Cu(111) and Ni(111). The minimal mismatch in lattice structure, under 5%, make these substrates ideal CVD candidates [32]. Graphene-metal interactions are often studied using density functional theory (DFT) calculations. From [32], which uses the local density approximation (LDA) method, we obtain values for the binding energy and interlayer distance, as shown in table 1. Graphene wrinkles observed with scanning tunnelling microscopy (STM) exhibit a wrinkle density of approximately 1 per μm 2 [15]. Subsequently, we select our half-length, L, as 500 nm.

Results-Model I: Arch wrinkles
In order to determine the geometric structures for arch wrinkles, we first need to solve (11) for the elliptic modulus k. As shown in figure 3, for a given value of l, there is a corresponding value of the elliptic modulus k that can be determined numerically.
Using the value of k, we can then plot the wrinkle profiles based on solutions (9) and (10). Figures 4 and 5 demonstrate the wrinkle structures for different nondimensional compression lengths l for graphene grown on Table 1. Interfacial properties between graphene and selected substrates obtained from DFT calculations [32]. Interaction between graphene layers is taken from [33].

Interlayer distance (Å)
Binding energy per unit area (eV Å −2 ) C-Ni 2.018 0.0913 C-Cu 3.260 0.0132 C-C 3.340 0.0214 substrates Ni and Cu, respectively. We can see from these figures that for the same value of compression length l, wrinkles are more pronounced for graphene grown Cu substrate than those of Ni substrate. This is because the binding energy for graphene and Ni is stronger than that of graphene and Cu. As such, wrinkles are more common for graphene grown on Cu. This may be due to the matching of hexagonal lattice of graphene and that of the Ni substrate leading to a stronger interaction between graphene and Ni [34]. The binding energies used in this study are shown in table 1. We also comment that for = 6 l , the conformations shown in figures 4 and 5 are not physical as the graphene-graphene separation distance at the neck of the wrinkle is less than the graphenegraphene equilibrium distance of 3.34Å. For such values of the non-dimensional compression length, Model II is applicable. However, they have been included ( = 6 l ) in the plots for mathematical interest. As shown in figures 6 and 7, we require more energy for the wrinkles to form. After they have formed, then less energy is required to compress the wrinkles further. For a given compression length, the energy stored in the wrinkle is higher on a Ni substrate than the same compression on a Cu substrate. As the compression length increases then at some point the C-C interaction from the left and right sides of the wrinkle becomes important and we must use Model II to capture all of the significant interactions. The compression length when the closest distance equals the graphene interlayer distance is given in table 2. As shown in this table the maximum compression length depends on the choice of substrate and the bending rigidity.

Results-Model II: Self-adhered wrinkles
For Model II we consider two separate curves  1 and  2 , of which only  2 depends on the choice of substrate. The conformation of the wrinkle section modelled by the curve  1 is driven by the C-C vdW interactions and thus is independent of the substrate. In figure 8 we show  1 for various values of the bending rigidity. However the curve  2 is dominated by graphene-substrate interaction. In figures 9 and 10 we show the conformation of the bottom of the self-adhered wrinkle  2 on a Ni and Cu substrate, respectively. The Ni substrate results in a tighter curve due to the stronger binding energy for graphene-Ni compared to graphene-Cu. Figures 11 and 12 show the relationship between the energy per unit length and the compression length for Ni and Cu substrates, respectively. The relative binding energy for graphene-Ni is stronger than graphenegraphene and we see in figure 11 that the energy increases with compression length. However the relative binding energy for graphene-Cu is weaker than graphene-graphene and consequently the energy decreases with compression length. As a result, the present model indicates that wrinkles would form more easily on a Cu substrate than a Ni substrate.
In table 3 we show the minimum compression length ℓ min for which Model II would apply. We note that this depends on the choice of substrate and bending rigidity. We further comment that there is a gap between the maximum compression length for Model I and the minimum compression length for Model II. This is due to a transitional conformation between Models I and II which has not been captured in the present work. Finally, we note that the bending rigidity of graphene does not depend on the growth substrate. It is the property of graphene which depends on temperature, shape/size of the sheet [35][36][37]. In particular, the bending rigidity of graphene decreases exponentially with increasing temperature [37], while increases with increasing size [35,36]. There is also a relation between ripples and the bending rigidity [35] and as further shown in this paper, the shape of wrinkle depends on the ratio of the bending rigidity of graphene and the van der Waals energy between graphene and substrate.

Conclusions
This paper investigates wrinkle structures of graphene formed on nickel and copper substrates during the CVD process. Due to the assumed symmetry of the problem, we model wrinkles as a curve in the xy-plane. Using the calculus of variations we derive a functional based on a competition of energy between the elastic energy stored in the curved graphene and the van der Waals interaction of the graphene sheet and a metal substrate. This model (Model I) is sufficient to examine graphene ripples and wrinkles before self-adhesion takes place. The model is then extended in Model II to account for the graphene-graphene van der Waals interactions in a selfadhered wrinkle. In both cases the variational problem can be solved analytically from which we may predict the geometry and energy profile for the resulting graphene conformation. In examining the interaction with different substrate we have chosen two cases: copper, where the graphenegraphene interaction is stronger than the graphene-copper interaction; and, nickel, where the graphenegraphene interaction is weaker than the graphene-nickel interaction. This disparity is the reason that the energy per unit length plot for Model II has a negative gradient in the copper case and a positive gradient in the nickel case. This issue also comes into play when comparing the energy per unit length plots between Model I and  Model II. With reference to figure 13 we see that for copper there is a clear point of intersection in the energy per unit length plots of Model I and Model II for all of the values of γ considered. This indicates that the self-adhered configuration becomes energetically favourable for a compression length which is valid for both models. This is in contrast to figure 14 which shows that for nickel, there are no points of intersection between the energy per unit length plots for Model I and Model II. This indicates that a third model is required that can be applied in the   domain of compression length between the maximum allowed by Model I and the minimum allowed by Model II. This presents a future direction for research. CVD grown wrinkles on copper substrate show height ranging from 20 Å to 60 Å examined using atomic force microscopy (AFM) [2,15,30,38], atomistic Monte Carlo simulations [39] and transmission electron microscopy (TEM) [40]. This correlates well with our theoretical data as shown in table 4. Particularly, Chaitoglou and Bertran [38] study the profile of the ripples formed in the graphene surface taken from their CVD experiments on copper substrate at temperature ranging from 970°C to 1070°C. They find the height of ripples to be 5.3±1nm at 1070°C. This agrees with our result of 54.97 Å (or 5.497 nm) shown in table 4 for copper substrate at the highest compression length. However, experimentally measured wrinkles are an order of magnitude wider than theoretical calculations [30]. It has been speculated that a collapsed self-adhered wrinkle could account for this discrepancy [2]. Regardless, this remains another issue requiring further theoretical investigation.
In conclusion, we have developed two models for arch-shaped and self-adhered wrinkles. In both cases, the principle parameter is the ratio between bending rigidity of graphene and the strength of the graphene-substrate van der Waals interaction. We find that less pronounced wrinkles are predicted for substrate with strong graphene-substrate van der Waals interaction. Furthermore, our model predicts the conformation of graphene wrinkle including key parameters, such as the curvature, width and height of the wrinkles. Our model also  predicts the transition between arch-shaped wrinkle and self-adhered wrinkle. Although, we perform our calculations on Ni and Cu solid substrates, our model is applicable to other solid metal substrates by employing appropriate values of the van der Waals strength. Future work includes catering for liquid metal/alloy substrates, trapped H 2 effect, and collapsed wrinkle conformation.