Over the past few decades, the measurement and analysis of thermal undulations has provided a route to estimate the mechanical properties of membranes. Theoretically, fluctuating elastic membranes have been studied mostly by Fourier analysis coupled with perturbation theory (to capture anharmonic effects), or by computer simulations of triangulated surfaces. These techniques as well as molecular dynamic simulations have also been used to study the thermal fluctuations of graphene. Here, we present a semi-analytic approach in which we view graphene as a triangulated membrane, but compute the statistical mechanical quantities using Gaussian integrals. The nonlinear coupling of in-plane strains with out-of-plane deflections is captured using a penalty energy. We recover well-known results for the scaling of the fluctuations with membrane size, but we show that the fluctuation profile strongly depends on boundary conditions and type of loading applied on the membrane. Our method quantitatively predicts the dependence of the thermal expansion coefficient of graphene on temperature and shows that it agrees with several experiments. We also make falsifiable predictions for the dependence of thermal expansion coefficient and the heat capacity of graphene on applied loads and temperature.

## Introduction

The thermal fluctuations of membranes, as well as sheets, plates, films, and surfaces, have been of interest for a few decades [1,2]. The problem has been studied experimentally [3,4], analytically [5–7], and through simulations [7,8]. A wealth of early literature on the subject focused on lipid bilayer membranes, polymer membranes, tethered surfaces, the crumpling transition, liquid crystals, etc. [3,9–11], which provide a ground for the interplay of geometry and statistical mechanics. In simulations, triangulated surface models combined with Monte Carlo and other simulation methods (such as time-dependent Ginzburg Landau) [12–16] have been the workhorse for the study of membrane statistical mechanics. Recently, we approached this problem using a different semi-analytic technique [17]. We discretized the membrane into equilateral triangle elements, just as has been done in earlier work, but confined ourselves to moderately large deflections so that the energy can be written as a quadratic form. We then used Gaussian integrals, instead of widely used simulation methods, such as Monte Carlo [12–14,18], to compute the partition function and free energy of this membrane. We applied our method to lipid bilayer membranes, which are liquid in-plane, and recovered classic results on the relation between tension and projected area that go back to Helfrich [5,9]. We also showed that our method can be applied to membranes with nonzero spontaneous curvature and heterogeneous bending modulus.

In this paper, the goal is to extend our semi-analytic method based on Gaussian integrals to solid membranes, or fluctuating elastic plates. The main difficulty in doing so is that the out-of-plane deflections of a solid membrane are coupled with the in-plane strains through the compatibility condition in von Karman plate theory [19]. Earlier simulation work on triangulated surface models as well as analytic theory has shown that undulations with nonzero Gauss curvature are suppressed because of this constraint [6,20–23]. In a similar vein, we show that a membrane explores fewer configurations if we enforce the compatibility constraint on our discretized plate through a penalty energy. The coupling between out-of-plane and in-plane displacement in solid membranes has also been analyzed in Fourier space using perturbation theory to get the anharmonic part of the spectrum [7,8,21]. To the best of our knowledge, none of these methods explore the effect of boundary conditions on the fluctuations. For example, analogous to filaments [24,25], we expect that a plate clamped on all its edges will suffer smaller fluctuations than one that is simply supported on all its edges. Our method allows us to explore these possibilities and consider fluctuations in the presence of shear loads. Note that the liquid membranes we studied before could not support shears. Also, our method can handle different membrane geometries. For example, in Ref. [17], we considered a square membrane with a circular patch in the middle that had a nonzero spontaneous curvature. Due to the nonzero spontaneous curvature, this patch assumed the shape of a spherical cap under zero loads. We analyzed the thermal fluctuations of this membrane and showed that our results were in agreement with earlier simulations using very different methods [26]. In these calculations, we did not need to enforce the compatibility constraint since the membrane was liquid in plane. To analyze the fluctuations of curved plates (or shells) using our method, it will be necessary to account for the compatibility constraint in curvilinear coordinates. While this is beyond the scope of our paper, progress has been made in the study of thermal fluctuations of shells [8] by triangulating them in a similar fashion as we do here.

As an illustration of the capability of our methods, we focus attention on the thermal fluctuations of graphene. The outstanding thermal, mechanical, and electrical properties of graphene [27,28] have received a lot of attention in the past decade, but its thermal fluctuations have received relatively less attention [21]. It is appreciated, however, that ripples on graphene have some bearing on its physical properties [21,29–32]. One such effect is its negative thermal expansion coefficient which has been studied through experiment and simulation [4,33,34]. Here, we quantitatively explain the dependence of the thermal expansion coefficient on temperature using our semi-analytic method. We also make predictions for the dependence of thermal expansion coefficient on loading (hydrostatic and shear) and temperature.

## Fluctuating Elastic Plate

### Von Karman Energy.

*xy*-plane in the reference (undeformed) configuration. We assume further that our plate is isotropic and use the von Karman plate energy expression [19], which consists of in-plane strain energy and bending energy as follows:

*Y*and

*ν*are the in-plane Young's modulus and Poisson's ratio of the material, respectively,

*K*is a mean curvature bending modulus, and

_{b}*K*is the Gauss curvature bending modulus.

_{G}*E*is the energy of in-plane strains, and

_{s}*E*is the energy of out-of-plane bending. The variables $\epsilon x,\u2009\epsilon y,\u2009\epsilon xy$, and

_{b}*w*in the integral are, respectively, the in-plane strains and out-plane displacement of the neutral plane. The subindices $\u2009,xx,,xy,,yy$ represent second derivatives, respectively, with respect to the reference coordinates

*x*and

*y*. We recall that the strains ${\epsilon x,\epsilon y,\epsilon xy}$ are not independent since

*F*along all the edges of our fluctuating plate. Then, from Eq. (2), the potential energy due to this tension is given by

*w*(

*x*,

*y*) (subject to the compatibility constraint equation (3)). It is instructive to make some estimates about the effect of fluctuations before we proceed further. From Eq. (1), we see that $Es\u223cYA\epsilon 2/2$, where

*ε*represents in-plane strains, and $Eb\u223cKbA\kappa 2/2$, where

*κ*represents a curvature. Using the equipartition theorem of statistical mechanics, we can estimate the mean-square fluctuations in these quantities as $\u2329\epsilon 2\u232a=kBT/YA$ and $\u2329\kappa 2\u232a=kBT/KbA$. If

*h*is the thickness of our plate, then the ratio of mean-square fluctuations in the bending strain to in-plane strain is $\u2329\epsilon 2\u232a/\u2329\kappa 2\u232ah2=Kb/Yh2$. For a graphene sheet,

*Y*≈ 1 TPa × 0.3 nm,

*K*≈ 10

_{b}^{−19}N·m, and

*h*≈ 0.3 nm [28], so that $Kb/Yh2\u224810\u22123$, and thus, in-plane strain fluctuations can be neglected in comparison to out-of-plane bending fluctuations. Due to the very high in-plane stretching modulus, we can also neglect the term $F(\epsilon x+\epsilon y)$ in Eq. (4) above. In the end, we assume that the plate has constant in-plane strains (no fluctuation) due to hydrostatic tension

*F*, which contributes a constant term

*C*to the expression for energy below

We will use this energy expression to compute the fluctuations of our membrane, remembering that *C* will have no effect on the fluctuations. To do so, we must calculate the partition function.

### Partition Function.

To compute the partition function, we will first discretize the plate and express the above energy as a function of node variables *w*. We discretize the plate into approximately $4N2/3$ equilateral triangles, as shown in Fig. 1(a), with $P\u22482N(N+1)/3$ node points, where *N* = *L*/*l*, in which *L* is the length of the square plate and *l* is the side of each triangle. Such a discretization has already been applied to fluid and solid membranes to study their statistical mechanics using Monte Carlo and other simulation techniques [37,38]. However, we compute the partition function using a different technique described in Ref. [17]. We approximate the discretized energy as quadratic expression in *w _{i}*, viz., $E=wMwT$, where the vector $w=[w1,w2,\u2026,wP]$ contains all the node displacements. The matrix

**M**is the stiffness matrix, and it is a function of $Kb,KG,F,L,\u2009and\u2009l$. The probability of finding the membrane in a given configuration of energy

*E*is proportional to $exp\u2009(\u2212E/kBT)$, where

*k*is the Boltzmann constant and

_{B}*T*is the absolute temperature. The partition function

*Z*is obtained by integration of $exp\u2009(\u2212E/kBT)$ over all allowed configurations of the plate. This integration is easy to carry out if the energy is a quadratic function of the displacements as shown in Refs. [25] and [39–42]. The partition function scales inversely with the square root of the determinant of matrix

**M**, or the square root of the product of all eigenvalues of

**M**as

### Enforcing the Compatibility Constraint.

*u*,

*v*, and

*w*are linear functions of the reference coordinates, and the compatibility equation is trivially satisfied inside every triangular element. However, at each node, the strains vary discontinuously and the Gaussian curvature is nonzero. Hence, in order to enforce the compatibility constraint at each node, we will penalize a violation of Eq. (3) by a large energy cost. For convenience, we define a compatibility function

*f*(

_{c}*x*,

*y*) by a scalar and integrating over an area around each node (the Voronoi cell). The total energy can be written as

*λ*(

*x*,

*y*) is a scalar with energy units, and

*A*is the area within the loop in Fig. 1(b). The summation runs over all nodes

_{i}*i*. In order to take advantage of Eq. (6), we take

*λ*to be a constant

*λ*within area

_{i}*A*such that it could come out of each integral. Then, we express $\u222cfcidAi$ as a quadratic expression in terms of the displacement vector

_{i}**w**just as we have done for the plate energy

*E*

**w**which consists of nodal displacement of node

*i*and six adjacent nodes as shown in Fig. 1(b). $Mic$ is a stiffness matrix of compatibility at node

*i*, such that $\u222cfcidAi=wiMicwiT$. $Mic$ will be computed in Sec. 2.4. Note that we can verify the effect of imposing compatibility at node

*i*by computing the average value of $wiMicwiT$. This is done by taking the partial derivative of the Gibb's free energy with respect to

*λ*as in Refs. [41] and [42]

_{i}### Compatibility Equation in Terms of Nodal Variables.

*γ*is the angle of the triangle element

_{i}*i*sharing this node. Let us assume that our plate is flat initially. We will now compute the angle of one element in terms of arbitrary nonzero node displacements

*w*,

*w*

_{1}, and

*w*

_{2}as shown in Fig. 1(b). Suppose the original length of the side of each element is

*l*. Then, the deformed length of the three sides is related to

*w*,

*w*

_{1}, and

*w*

_{2}as

*w*and get the integral of the Gaussian curvature at a node as

*f*, could be evaluated up to quadratic order in the displacements $wi$ at the nodes, governed by matrix $Mic$ in Eq. (9). We will now choose all

_{c}*λ*equal to a same value

_{i}*λ*because there is no reason to penalize two nodes differently. Then, the energy expression is

*λ*contributes to the fluctuation through the boundaries. Note that in Fig. 1(b), the integration contours around two adjacent nodes share one edge, which is traversed in opposite directions. Therefore, if we choose

*λ*to be the same for these two nodes, then the integration in Eq. (18) along this edge will cancel due to opposite signs of the contribution from the two contours that share it. Similarly, if

*λ*is chosen to be the same for all nodes, this cancelation will happen for all edges in the interior of the plate, and the entire energy penalty due to the compatibility constraint goes to the boundary of the plate. This is not surprising if one remembers the Gauss–Bonnet theorem for the integral of the Gaussian curvature [44]. Thus, the summation of the energy cost of violating the compatibility constraint at each node can be written as a quadratic form

where the matrix $Mc$ has nonzero entries only for the boundary nodes. Furthermore, we find that $Mc$ is positive semidefinite irrespective of the boundary conditions, so that we only need to consider *λ* > 0 in Eq. (19). Now, in order to enforce the compatibility constraint, we choose *λ* sufficiently large so that the variance of out-of-plane deflections *w*(*x*, *y*) is independent of *λ*. In reality, by doing so, we have required that *f _{c}* is close to zero on average over the whole plate. In order to justify this weaker constraint, we will compute $\u2329\u222cfcidAi\u232a$ over a Voronoi cell surrounding each node to make sure that they are approaching zero as well. This is exactly the quantity in Eq. (10). We examine it at every node

*i*for a solid membrane and a fluid membrane in Fig. 2 for two different boundary conditions as given in Sec. 3.1. Since this quantity is very small (on the order of 10

^{−6}) at all nodes for the solid membrane, we have ensured that the compatibility equation is satisfied everywhere. In contrast, for a liquid membrane with the same mechanical properties and boundary conditions,

*f*is much larger as shown in Fig. 2. We also show (a) how the variance of the fluctuations of a solid membrane in Figs. 3(a) and 3(b) becomes independent of

_{c}*λ*as it becomes sufficiently large (see the solid lines associated with left axis), and (b) how the average value of

*f*over all nodes approaches zero as

_{c}*λ*becomes sufficiently large (see the dashed line associated with right axis). We see that

*λ*= 10

^{9}pN nm is a large enough value and use it in all our subsequent calculations of the variance. The stage is now set to examine the effect of boundary conditions on a fluctuating graphene sheet.

## Results

### Analysis of Variance and Effect of Boundary Conditions.

From this, we can get the variance of the out-plane deflection $\u2329wi2\u232a$ as a function of reference position on our membrane. We choose membrane size, bending moduli, tension, boundary condition, and finite element discretization as given in Table 1. First, we examine the overall variance profile of the whole plate with two different boundary conditions—(a) one edge hinged, three edges free and (b) two opposite edges hinged, two edges free. The variance is plotted as function of node position in the reference configuration shown in Figs. 3(a) and 3(b). The results correspond to computation groups A and B in Table 1. In order to show the effect of the compatibility constraint, we plot the variance profile of a fluid membrane with the same *K _{b}* and

*K*in each plot for comparison. Recall that since a fluid membrane has zero in-plane shear modulus, we can choose $\epsilon xy(x,y)$ to satisfy the compatibility equation for a given out-of-plane deflection profile

_{G}*w*(

*x*,

*y*) and in-plane strain profile $\epsilon x(x,y),\u2009\epsilon y(x,y)$. The compatibility constraint reduces the thermal fluctuations of solid membranes compared to solid membranes for the same boundary conditions.

Group | L or L_{1} (nm) | K (pN nm)_{b} | K (pN nm)_{G} | F (pN nm^{−1}) | Boundary condition | N |
---|---|---|---|---|---|---|

A | 1000 | 82 | −53.3 | 0.001 | One end hinged | 50 |

B | 1000 | 82 | −53.3 | 0.001 | Two ends hinged | 50 |

C | 5–150 | 41 | −26.65 | 0 | One end hinged | 50 |

D | 5–150 | 41 | −26.65 | 0 | Two ends hinged | 50 |

E | 5–150 | 82 | −53.3 | 0 | One end hinged | 50 |

F | 5–150 | 82 | −53.3 | 0 | Two ends hinged | 50 |

G | 5–150 | 205 | −133.25 | 0 | One end hinged | 50 |

H | 5–150 | 205 | −133.25 | 0 | Two ends hinged | 50 |

I | 5–150 | 3.28 | −2.132 | 0 | One end hinged | 50 |

J | 1000 (= L_{2}) | 82 | −53.3 | 0 | Shear S = 0.001 pN nm^{−1} | 50 |

K | 1000 (= L_{2}) | 82 | −53.3 | 0 | Shear S = 0.01 pN nm^{−1} | 50 |

L | 1000 (L_{2} = 500) | 28.7 | −18.655 | 0 | Shear S = 0.01 pN nm^{−1} | 50 |

Group | L or L_{1} (nm) | K (pN nm)_{b} | K (pN nm)_{G} | F (pN nm^{−1}) | Boundary condition | N |
---|---|---|---|---|---|---|

A | 1000 | 82 | −53.3 | 0.001 | One end hinged | 50 |

B | 1000 | 82 | −53.3 | 0.001 | Two ends hinged | 50 |

C | 5–150 | 41 | −26.65 | 0 | One end hinged | 50 |

D | 5–150 | 41 | −26.65 | 0 | Two ends hinged | 50 |

E | 5–150 | 82 | −53.3 | 0 | One end hinged | 50 |

F | 5–150 | 82 | −53.3 | 0 | Two ends hinged | 50 |

G | 5–150 | 205 | −133.25 | 0 | One end hinged | 50 |

H | 5–150 | 205 | −133.25 | 0 | Two ends hinged | 50 |

I | 5–150 | 3.28 | −2.132 | 0 | One end hinged | 50 |

J | 1000 (= L_{2}) | 82 | −53.3 | 0 | Shear S = 0.001 pN nm^{−1} | 50 |

K | 1000 (= L_{2}) | 82 | −53.3 | 0 | Shear S = 0.01 pN nm^{−1} | 50 |

L | 1000 (L_{2} = 500) | 28.7 | −18.655 | 0 | Shear S = 0.01 pN nm^{−1} | 50 |

Second, we examine the relation between variance and membrane size while fixing the mechanical properties. An earlier theoretical result gives that the average fluctuation $w\xaf$ (square root of variance) scales with some power of membrane size *L*, or as *L ^{δ}*.

*δ*is proved to be 1 for a fluid membranes as in Refs. [6] and [21] and around 0.6–0.8 as given in Refs. [6] and [21–23] for solid membranes. We recover this power-law scaling in our calculations as shown in Fig. 3(d). The parameters of our membrane are chosen from groups C to H in Table 1. However, we find that the power

*δ*for a given solid membrane depends on the boundary conditions. When only one edge is hinged and the others are free,

*δ*= 0.8321, and when two opposite edges are hinged and two are free,

*δ*= 0.5968. These powers are close to the range given in Refs. [6] and [21].

The fluctuating plate model studied in our paper is neither this ideal membrane (with *K _{b}* = 0 and

*μ*≠ 0) nor a fluid membrane (with

*K*≠ 0 and

_{b}*μ*= 0). Therefore, the resulting power is between 0.5 and 1, and it depends on the applied boundary condition. In order to further verify our model, we give results from another calculation in the inset of Fig. 3(d) using parameters from group I in Table 1. We cannot choose the bending modulus

*K*and Gaussian curvature modulus

_{b}*K*to be exactly zero since this results in a singular stiffness matrix

_{G}**M**. So, we choose

*K*and

_{b}*K*to be smaller than the thermal energy scale 1

_{G}*k*. We find in Fig. 3(d) inset that the power

_{B}T*δ*for this membrane is

*δ*= 0.5542, which is closer to the theory in Refs. [5] and [45].

### Study of Ripples Under Shear Loading.

*F*with some other boundary condition, for example, a shear loading

*S*as in Ref. [46]. Along one pair of opposite edges, we apply a tension of magnitude +

*S*, while along the other pair of opposite edges, we apply a compression −

*S*. This will result in shear loading on the plate. To account for this loading, Eq. (4) is changed to

*w*(

*x*,

*y*) of our plate can be expressed as a double Fourier series [5]

*q*is the amplitude of the normal mode corresponding to

_{mn}*m*,

*n*. Note that we assume our plate to be rectangular with sides

*L*

_{1}in

*x*direction and

*L*

_{2}in

*y*direction. When we plug the above expression into Eq. (24), we get

*m*,

*n*) should possess energy

*k*/2. Therefore, we can compute the amplitude of each mode as

_{B}T*m*,

*n*) with maximum amplitude should be given by

where [⋅] gives the nearest positive integer. The mode shapes observed in Figs. 4(a)–4(c) are consistent with the result of plugging parameters of groups J–L into Eq. (28). The mode making the dominant contribution to the variance changes with the loading because the solid membranes can support shear. In contrast, in fluid membranes, the dominant contribution to the variance always comes from the lowest normal mode *n* = 1, *m* = 1 because they can only sustain hydrostatic tension. It is trivial to extend these calculations to combined hydrostatic and shear loading on the plate. It is evident from the above results that different loading and boundary conditions will result in different ripple profiles in fluctuating solid plates.

### Quantitative Analysis of the Negative Thermal Expansion Coefficient of Graphene.

^{−6}K

^{−1}. This remarkable feature of graphene has been measured and predicted in Refs. [4,21,33,34], and [48]. However, all these works were based either on experiment or molecular simulation. Here, we give a quantitative explanation for the negative thermal expansion coefficient using our fluctuating elastic plate model. The Gibbs free energy

*G*(

*F*,

*T*) of the whole plate is related to the partition function

*Z*as $G=\u2212kBTlnZ$. We can compute the properties of the plate by computing derivatives of the free energy. First, the reduction in projected area caused by thermal fluctuation is given by $\Delta A=A(\u221e,T)\u2212A(F,T)$, where

*A*(

*∞*,

*T*) is the membrane area at very large tension such that all the undulations are stretched out and it is flat. The area reduction is conjugate to the tension

*F*in the free energy of the membrane and can be computed as Δ

*A*= −

*∂G*/

*∂F*. Since all the terms related to the force

*F*in the expression for the partition function

*Z*are included in the matrix

**M**, we have

*A*is the original area of the membrane at zero tension and zero temperature. The increase in out-of-plane bending deflections of the plate at higher temperatures reduces the projected area and results in a negative thermal expansion coefficient. The dependence of the bending modulus

*K*on the temperature also contributes to the thermal expansion coefficient. This dependence has been studied in Ref. [21], and we summarize it through the following fit:

_{b}*K*and

_{b}*T*are pN nm and K, respectively. We evaluate Δ

*A*/

*A*at

*F*= 0.001 pN/nm of a 100 nm × 100 nm graphene sheet, for various values of

*T*. By differentiating Eq. (29) with respect to temperature

*T*, we get the thermal expansion coefficient

We plot our result in Fig. 5(a) together with other simulation and experimental results. It is apparent that our method agrees quite well with the earlier work. In particular, we find that the thermal expansion coefficient increases with increasing temperature as has been found in several experiments [4,33,34]. Although the experimental data for *α* is available over a range 200–400 K, we have shown in the inset of Fig. 5(a) that *α* < 0 over a much broader range up to 800 K. We also must not lose sight of the fact that this figure assumes that the physics behind the negative *α* is bending fluctuations, which may not be the right physics for a solid near 0 K. We also carefully examine the effect of distinct hydrostatic tension and shear loading on thermal expansion coefficient as given in Fig. 5(b). We see that hydrostatic tension does not affect the thermal expansion coefficient much. The shear loading, however, enormously affects it. First, as the shear load increases, the slope with respect to temperature *T* changes from positive to negative. Second, at shear loading equal to 0.001 pN nm^{−1}, the thermal expansion coefficient becomes positive. Note also that the slope of the *α* versus *T* curve is related to the heat capacity *C _{F}* = −

*T∂*

^{2}

*G*/

*∂T*

^{2}at given tension

*F*or shear loading

*S*through

## Conclusion

In this paper, we demonstrate a newly developed semi-analytic method to study the thermal fluctuations of solid membranes. We discretized the membranes using triangular elements and write their energy as a quadratic form starting from von Karman plate theory. We then evaluate the partition function using Gaussian integrals. We account for the nonlinear coupling of in-plane strains and out-of-plane deflections using a penalty method. Once the partition function of the fluctuating membrane is known, several thermodynamic quantities can be determined by calculating its derivatives. We have utilized this idea to illuminate how loading and boundary conditions affect the fluctuation profile (or ripples) of membranes. We have also shown that our method can quantitatively explain the dependence of the negative thermal expansion coefficient of graphene on temperature. We have made predictions for how shear loads can affect the thermal expansion coefficient that can be tested in experiments. Our method can be used to make quantitative predictions for other two-dimensional materials.

## Acknowledgment

We acknowledge support for this work through an NSF CAREER Award Grant No. NSF CMMI 0953548.