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 [57], 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,911], 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) [1216] 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 [1214,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,2023]. 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,2932]. 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.

In order to study the thermal fluctuations of a single sheet of graphene, we must first write an appropriate energy expression. Assuming that number of atoms comprising the sheet is large, we will approximate it as an elastic plate capable of moderately large deflections [28,35,36]. The plate is assumed to lie on the 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:
$Evon=Es+Eb=∬dxdyY2[(εx+εy)2−2(1−ν)(εxεy−εxy2)]+∬dxdy[Kb2(w,xx+w,yy)2+KG(w,xxw,yy−w,xy2)]$
(1)
where Y and ν are the in-plane Young's modulus and Poisson's ratio of the material, respectively, Kb is a mean curvature bending modulus, and KG is the Gauss curvature bending modulus. Es is the energy of in-plane strains, and Eb is the energy of out-of-plane bending. The variables $εx, εy, εxy$, and w in the integral are, respectively, the in-plane strains and out-plane displacement of the neutral plane. The subindices $,xx,,xy,,yy$ represent second derivatives, respectively, with respect to the reference coordinates x and y. We recall that the strains ${εx,εy,εxy}$ are not independent since
$εx=u,x+w,x22εy=v,y+w,y22εxy=u,y+v,x2+w,xw,y2$
(2)
These strain–displacement relations are connected through the compatibility equation
$εx,yy+εy,xx−2εxy,xy=w,xy2−w,xxw,yy$
(3)
Now, consider applying a two-dimensional hydrostatic tension of magnitude F along all the edges of our fluctuating plate. Then, from Eq. (2), the potential energy due to this tension is given by
$Ef=−∬F(u,x+v,y)dxdy=−∬F(εx−w,x22+εy−w,y22)dxdy$
(4)
Equations (1) and (4) together give the energy of a configuration of the plate which is determined by $εx(x,y), εy(x,y), εxy(x,y)$, and 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∼YAε2/2$, where ε represents in-plane strains, and $Eb∼KbAκ2/2$, where κ represents a curvature. Using the equipartition theorem of statistical mechanics, we can estimate the mean-square fluctuations in these quantities as $〈ε2〉=kBT/YA$ and $〈κ2〉=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 $〈ε2〉/〈κ2〉h2=Kb/Yh2$. For a graphene sheet, Y ≈ 1 TPa × 0.3 nm, Kb ≈ 10−19 N·m, and h ≈ 0.3 nm [28], so that $Kb/Yh2≈10−3$, 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(εx+ε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
$E=Eb+Ef+C=∬dxdy[Kb2(w,xx+w,yy)2+KG(w,xxw,yy−w,xy2)+F(w,x22+w,y22)]+C$
(5)

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≈2N(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 wi, viz., $E=wMwT$, where the vector $w=[w1,w2,…,wP]$ contains all the node displacements. The matrix M is the stiffness matrix, and it is a function of $Kb,KG,F,L, and l$. The probability of finding the membrane in a given configuration of energy E is proportional to $exp (−E/kBT)$, where kB is the Boltzmann constant and T is the absolute temperature. The partition function Z is obtained by integration of $exp (−E/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 [3942]. 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

$Z=(2πkBT)PdetM$
(6)

### Enforcing the Compatibility Constraint.

In order to extend this technique to solid plates, we must carry out the partition sum only over those configurations that satisfy the compatibility equation (3). This equation constrains the displacements at every point on the plate. Now, the question is how do we enforce this constraint everywhere on the plate. Recall that we discretize our plate into equilateral triangle elements in the reference or undeformed configuration. When the plate is deformed, the elements remain triangular and flat with straight edges. Hence, all the displacements 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
$fc=εx,yy+εy,xx−2εxy,xy−w,xy2+w,xxw,yy$
(7)
The energy cost can be computed by multiplying fc(x, y) by a scalar and integrating over an area around each node (the Voronoi cell). The total energy can be written as
$E¯=E+∬λfcdA=E+∑i∬λifcidAi$
(8)
where λ(x, y) is a scalar with energy units, and Ai is the area within the loop in Fig. 1(b). The summation runs over all nodes i. In order to take advantage of Eq. (6), we take λ to be a constant λi within area Ai such that it could come out of each integral. Then, we express $∬fcidAi$ as a quadratic expression in terms of the displacement vector w just as we have done for the plate energy E
$E¯=wMwT+∑iλiwiMicwiT$
(9)
Here, $wi=[w,w1,w2,w3,w4,w5,w6]$ is a subset of 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 $∬fcidAi=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 λi as in Refs. [41] and [42]
$〈wiMicwiT〉=∂G∂λi$
(10)

### Compatibility Equation in Terms of Nodal Variables.

A convenient way to calculate the integral of Gaussian curvature in RHS of Eq. (3) at a node is given by Magid et al. [43] as
$∬CGdA=2π−∑γi$
(11)
where γi is the angle of the triangle element 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, w1, and w2 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, w1, and w2 as
$a=(w−w1)2+l2b=(w−w2)2+l2c=(w1−w2)2+l2$
(12)
So, the deformed angle is given by
$γ=arccosa2+b2−c22ab$
(13)
We plug in Eq. (12) and then calculate the change of angle with respect to the original equilateral triangle. In order to carry out the integration for the partition function, we approximate this expression quadratically
$γ−π3=2(w1−w2)2−(w−w1)2−(w−w2)223l2$
(14)
We apply this calculation to other triangle elements sharing the node with out-of-plane displacement w and get the integral of the Gaussian curvature at a node as
$∬CGdA=∑i=16(π3−γi)=∑i=16(w−wi)2−∑i=16(wi−wi+1)23l2$
(15)
This takes care of the right-hand side of the compatibility constraint equation (3). The contribution to the elastic bending energy due to the Gaussian curvature can also be calculated using the above formula. Now, we focus on the left-hand side of Eq. (3), which involves second derivatives of the in-plane strains. To evaluate the integral of this quantity over a Voronoi cell surrounding a node, we choose the boundary of one such cell, which is a closed integration contour as shown in Fig. 1(b). Then, from Stokes theorem, we get
$I=∬(εx,yy+εy,xx−2εxy,xy)dA=∮(εxy,x−εx,y)dx+∮(εy,x−εxy,y)dy$
(16)
The terms $εxy,x$ and $εxy,y$ in the integrands will contribute zero because our integration contour is closed. Hence, we are left with
$I=∮εy,xdy−∮εx,ydx$
(17)
Now, we plug Eq. (2) into the above integral and similarly drop those terms in the integrands whose contribution will be zero due to the closed contour. Finally, we are left with
$I=∮[12∂∂x(∂w∂y)2dy−12∂∂y(∂w∂x)2dx]$
(18)
The expression above can be easily evaluated in terms of the nodal variables for each element. Thus, both sides of Eq. (3), or the compatibility function fc, 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 λi equal to a same value λ because there is no reason to penalize two nodes differently. Then, the energy expression is
$E¯=wMwT+λ∑iwiMicwiT$
(19)
We will now show that λ 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
$∑iwiMicwiT=wMcwT$
(20)

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 fc is close to zero on average over the whole plate. In order to justify this weaker constraint, we will compute $〈∬fcidAi〉$ 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, fc 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 λ as it becomes sufficiently large (see the solid lines associated with left axis), and (b) how the average value of fc over all nodes approaches zero as λ becomes sufficiently large (see the dashed line associated with right axis). We see that λ = 109 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.

One important result which follows from the calculation of the partition function is the covariance matrix $〈wiwj〉$ [25,4042], which can be calculated from the inverse of the stiffness matrix M as
$〈wiwj〉=kBT[M−1]ij$
(21)

From this, we can get the variance of the out-plane deflection $〈wi2〉$ 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 Kb and KG in each plot for comparison. Recall that since a fluid membrane has zero in-plane shear modulus, we can choose $εxy(x,y)$ to satisfy the compatibility equation for a given out-of-plane deflection profile w(x, y) and in-plane strain profile $εx(x,y), εy(x,y)$. The compatibility constraint reduces the thermal fluctuations of solid membranes compared to solid membranes for the same boundary conditions.

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¯$ (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 [2123] 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].

In the analysis given by Boal et al. [5,45], an ideal membrane with nonzero in-plane shear resistance μ, but zero bending resistance Kb, has been shown to have an average thermal fluctuation $w¯∝L0.5$ given by
$w¯=0.62L0.5(kBTY)0.25$
(22)

The fluctuating plate model studied in our paper is neither this ideal membrane (with Kb = 0 and μ ≠ 0) nor a fluid membrane (with Kb ≠ 0 and μ = 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 Kb and Gaussian curvature modulus KG to be exactly zero since this results in a singular stiffness matrix M. So, we choose Kb and KG to be smaller than the thermal energy scale 1 kBT. We find in Fig. 3(d) inset that the power δ for this membrane is δ = 0.5542, which is closer to the theory in Refs. [5] and [45].

An interesting possibility for us is to replace the hydrostatic tension 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
$Ef=−∬S(u,x−v,y)dxdy=−∬S(εx−w,x22−εy+w,y22)dxdy$
(23)
This results in a sign change in Eq. (5) as
$E=∬dxdy[Kb2(w,xx+w,yy)2+KG(w,xxw,yy−w,xy2)+S(w,x22−w,y22)]$
(24)
In Table 1, groups J–L, we give parameters for this type of problem. The resulting variance profiles for the out-of-plane deflection are given in Fig. 4 for three different combinations of parameters. All three membranes are simply supported on the four edges. Each combination of properties results in a distinct variance profile. Next, we are going to explain the relation between the parameters and the variance profile obtained from our computation. Due to hinged boundary conditions on all edges, the deflections w(x, y) of our plate can be expressed as a double Fourier series [5]

$w(x,y)=∑(m,n)qmn sin(nπxL1)sin(mπyL2)$
(25)
where qmn is the amplitude of the normal mode corresponding to m, n. Note that we assume our plate to be rectangular with sides L1 in x direction and L2 in y direction. When we plug the above expression into Eq. (24), we get
$E=∑(m,n)qmn2[Kbπ48L1L2(n2+m2)2+Sπ28(n2−m2)]$
(26)
From the equipartition theorem of statistical mechanics [47], each Fourier mode (m, n) should possess energy kBT/2. Therefore, we can compute the amplitude of each mode as
$〈qmn2〉=4kBTL1L2Kbπ4(n2+m2)2+SL1L2π2(n2−m2)$
(27)
Therefore, the mode (m, n) with maximum amplitude should be given by
$n=1; m=[1+2SL1L2Kbπ2−12]$
(28)

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.

We will use our model to analyze the negative thermal expansion coefficient of graphene, which is around −10−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=−kBTlnZ$. 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 $ΔA=A(∞,T)−A(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
$ΔAA=−kBT2A∂∂Fln(detM)$
(29)
where 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 Kb 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:
$Kb=131.36+200· tan hT1500$
(30)
where the units of Kb and 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
$α=−∂2∂T∂F[kBT2AlndetM]$
(31)

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 CF = −T∂2G/∂T2 at given tension F or shear loading S through

$∂α∂T=−1A∂3G∂F∂T2=1AT∂CF∂F$
(32)

In Fig. 5(b) top panel, we see that this slope changes sign at a certain value of S. The equation above shows that a measurement of the heat capacity as a function of S can be a good way of testing our predictions in addition to experiments of the type carried out in Refs. [4,33], and [34].

## 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.

## References

References
1.
Helfrich
,
W.
,
1973
, “
Elastic Properties of Lipid Bilayers: Theory and Possible Experiments
,”
Z. Naturforsch. C
,
28
(
11
), pp.
693
703
.
2.
Weeks
,
J. D.
,
1977
, “
Structure and Thermodynamics of the Liquid–Vapor Interface
,”
J. Chem. Phys.
,
67
(
7
), pp.
3106
3121
.
3.
Evans
,
E.
, and
Rawicz
,
W.
,
1990
, “
Entropy-Driven Tension and Bending Elasticity in Condensed-Fluid Membranes
,”
Phys. Rev. Lett.
,
64
(
17
), p.
2094
.
4.
Bao
,
W.
,
Miao
,
F.
,
Chen
,
Z.
,
Zhang
,
H.
,
Jang
,
W.
,
Dames
,
C.
, and
Lau
,
C. N.
,
2009
, “
Controlled Ripple Texturing of Suspended Graphene and Ultrathin Graphite Membranes
,”
Nat. Nanotechnol.
,
4
(
9
), pp.
562
566
.
5.
Boal
,
D.
,
2002
,
Mechanics of the Cell
,
Cambridge University Press
,
Cambridge, UK
.
6.
David
,
F.
,
Nelson
,
D.
,
Piran
,
T.
, and
Weinberg
,
S.
,
1989
,
Statistical Mechanics of Membranes and Surfaces
,
D.
Nelson
,
T.
Piran
, and
S.
Weinberg
, eds., World Scientific, River Edge, NJ.
7.
Kosmrlj
,
A.
, and
Nelson
,
D. R.
,
2014
, “
Thermal Excitations of Warped Membranes
,”
Phys. Rev. E
,
89
(
2
), p.
022126
.
8.
Paulose
,
J.
,
Vliegenthart
,
G. A.
,
Gompper
,
G.
, and
Nelson
,
D. R.
,
2012
, “
Fluctuating Shells Under Pressure
,”
,
109
(
48
), pp.
19551
19556
.
9.
Helfrich
,
W.
,
1975
, “
Out-of-Plane Fluctuations of Lipid Bilayers
,”
Z. Naturforsch. Sect. C: Biosci.
,
30
(
6
), p.
841
.
10.
Milner
,
S. T.
, and
Safran
,
S. A.
,
1987
, “
Dynamical Fluctuations of Droplet Microemulsions and Vesicles
,”
Phys. Rev. A
,
36
(
9
), p.
4371
.
11.
Rodriguez-Garcia
,
R.
,
Mell
,
M.
,
Lopez-Montero
,
I.
,
Netzel
,
J.
,
Hellweg
,
T.
, and
Monroy
,
F.
,
2011
, “
Polymersomes: Smart Vesicles of Tunable Rigidity and Permeability
,”
Soft Matter
,
7
(
4
), pp.
1532
1542
.
12.
Auth
,
T.
, and
Gompper
,
G.
,
2013
, “
Fluctuation Pressure of Biomembranes in Planar Confinement
,”
Phys. Rev. E
,
88
(
1
), p.
010701
.
13.
Ramakrishnan
,
N.
,
Kumar
,
P. S.
, and
Ipsen
,
J. H.
,
2010
, “
Monte Carlo Simulations of Fluid Vesicles With In-Plane Orientational Ordering
,”
Phys. Rev. E
,
81
(
4
), p.
041922
.
14.
Zhang
,
T.
,
Li
,
X.
, and
Gao
,
H.
,
2014
, “
,”
J. Mech. Phys. Solids
,
67
, pp.
2
13
.
15.
Harmandaris
,
V. A.
, and
Deserno
,
M.
,
2006
, “
A Novel Method for Measuring the Bending Rigidity of Model Lipid Membranes by Simulating Tethers
,”
J. Chem. Phys.
,
125
(
20
), p.
204905
.
16.
Ramakrishnan
,
N.
,
Kumar
,
P. S.
, and
,
R.
,
2014
, “
Mesoscale Computational Studies of Membrane Bilayer Remodeling by Curvature-Inducing Proteins
,”
Phys. Rep.
,
543
(
1
), pp.
1
60
.
17.
Liang
,
X.
, and
Purohit
,
K. P.
,
2015
, “
A Fluctuating Elastic Plate and a Cell Model for Lipid Membranes
,”
J. Mech. Phys. Solids
,
90
, pp.
29
44
.
18.
Hanlumyuang
,
Y.
,
Liu
,
L.
, and
Sharma
,
P.
,
2014
, “
Revisiting the Entropic Force Between Fluctuating Biological Membranes
,”
J. Mech. Phys. Solids
,
63
, pp.
179
186
.
19.
Audoly
,
B.
, and
Pomeau
,
Y.
,
2010
,
Elasticity and Geometry
,
Oxford University Press
,
Oxford, UK
.
20.
Nelson
,
D. R.
, and
Peliti
,
L.
,
1987
, “
Fluctuations in Membranes With Crystalline and Hexatic Order
,”
J. Phys.
,
48
(
7
), pp.
1085
1092
.
21.
Fasolino
,
A.
,
Los
,
J. H.
, and
Katsnelson
,
M. I.
,
2007
, “
Intrinsic Ripples in Graphene
,”
Nat. Mater.
,
6
(
11
), pp.
858
861
.
22.
Le Doussal
,
P.
, and
,
L.
,
1992
, “
Self-Consistent Theory of Polymerized Membranes
,”
Phys. Rev. Lett.
,
69
(
8
), p.
1209
.
23.
Abraham
,
F. F.
, and
Nelson
,
D. R.
,
1990
, “
Diffraction From Polymerized Membranes
,”
Science
,
249
(
4967
), pp.
393
397
.
24.
Purohit
,
P. K.
,
Arsenault
,
M. E.
,
Goldman
,
Y.
, and
Bau
,
H. H.
,
2008
, “
The Mechanics of Short Rod-Like Molecules in Tension
,”
Int. J. Non-Linear Mech.
,
43
(
10
), pp.
1056
1063
.
25.
Su
,
T.
, and
Purohit
,
P. K.
,
2010
, “
Thermomechanics of a Heterogeneous Fluctuating Chain
,”
J. Mech. Phys. Solids
,
58
(
2
), pp.
164
186
.
26.
Agrawal
,
N. J.
, and
,
R.
,
2009
, “
Calculation of Free Energies in Fluid Membranes Subject to Heterogeneous Curvature Fields
,”
Phys. Rev. E
,
80
(
1
), p.
011925
.
27.
Bunch
,
J. S.
,
Van Der Zande
,
A. M.
,
Verbridge
,
S. S.
,
Frank
,
I. W.
,
Tanenbaum
,
D. M.
,
Parpia
,
J. M.
, and
McEuen
,
P. L.
,
2007
, “
Electromechanical Resonators From Graphene Sheets
,”
Science
,
315
(
5811
), pp.
490
493
.
28.
Lee
,
C.
,
Wei
,
X.
,
Kysar
,
J. W.
, and
Hone
,
J.
,
2008
, “
Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene
,”
Science
,
321
(
5887
), pp.
385
388
.
29.
Los
,
J. H.
,
Katsnelson
,
M. I.
,
Yazyev
,
O. V.
,
Zakharchenko
,
K. V.
, and
Fasolino
,
A.
,
2009
, “
Scaling Properties of Flexible Membranes From Atomistic Simulations: Application to Graphene
,”
Phys. Rev. B
,
80
(
12
), p.
121405
.
30.
Zakharchenko
,
K. V.
,
Los
,
J. H.
,
Katsnelson
,
M. I.
, and
Fasolino
,
A.
,
2010
, “
Atomistic Simulations of Structural and Thermodynamic Properties of Bilayer Graphene
,”
Phys. Rev. B
,
81
(
23
), p.
235439
.
31.
Garcia-Sanchez
,
D.
,
van der Zande
,
A. M.
,
Paulo
,
A. S.
,
Lassagne
,
B.
,
McEuen
,
P. L.
, and
Bachtold
,
A.
,
2008
, “
Imaging Mechanical Vibrations in Suspended Graphene Sheets
,”
Nano Lett.
,
8
(
5
), pp.
1399
1403
.
32.
He
,
Y. Z.
,
Li
,
H.
,
Si
,
P. C.
,
Li
,
Y. F.
,
Yu
,
H. Q.
,
Zhang
,
X. Q.
, and
Liu
,
X. F.
,
2011
, “
Dynamic Ripples in Single Layer Graphene
,”
Appl. Phys. Lett.
,
98
(
6
), p.
063101
.
33.
Yoon
,
D.
,
Son
,
Y. W.
, and
Cheong
,
H.
,
2011
, “
Negative Thermal Expansion Coefficient of Graphene Measured by Raman Spectroscopy
,”
Nano Lett.
,
11
(
8
), pp.
3227
3231
.
34.
Pan
,
W.
,
Xiao
,
J.
,
Zhu
,
J.
,
Yu
,
C.
,
Zhang
,
G.
,
Ni
,
Z.
,
Watanabe
,
K.
,
Taniguchi
,
T.
,
Shi
,
Y.
, and
Wang
,
X.
,
2012
, “
Biaxial Compressive Strain Engineering in Graphene/Boron Nitride Heterostructures
,”
Sci. Rep.
,
2
, p. 893.
35.
Wei
,
X.
,
Fragneaud
,
B.
,
Marianetti
,
C. A.
, and
Kysar
,
J. W.
,
2009
, “
Nonlinear Elastic Behavior of Graphene: Ab Initio Calculations to Continuum Description
,”
Phys. Rev. B
,
80
(
20
), p.
205407
.
36.
Lee
,
G. H.
,
Cooper
,
R. C.
,
An
,
S. J.
,
Lee
,
S.
,
van der Zande
,
A.
,
Petrone
,
N.
, and
Kysar
,
J. W.
,
2013
, “
High-Strength Chemical-Vapor-Deposited Graphene and Grain Boundaries
,”
Science
,
340
(
6136
), pp.
1073
1076
.
37.
Gompper
,
G.
, and
Kroll
,
D. M.
,
1996
, “
Random Surface Discretizations and the Renormalization of the Bending Rigidity
,”
J. Phys. I
,
6
(
10
), pp.
1305
1320
.
38.
Fraternali
,
F.
, and
Marcelli
,
G.
,
2012
, “
A Multiscale Approach to the Elastic Moduli of Biomembrane Networks
,”
Biomech. Model. Mechanobiol.
,
11
(
7
), pp.
1097
1108
.
39.
Flory
,
P. J.
, and
Volkenstein
,
M.
,
1969
, “
Statistical Mechanics of Chain Molecules
,”
Biopolymers
,
8
(
5
), pp.
699
700
.
40.
Zhang
,
Y.
, and
Crothers
,
D. M.
,
2003
, “
Statistical Mechanics of Sequence-Dependent Circular DNA and Its Application for DNA Cyclization
,”
Biophys. J.
,
84
(
1
), pp.
136
153
.
41.
Su
,
T.
, and
Purohit
,
P. K.
,
2011
, “
Fluctuating Elastic Filaments Under Distributed Loads
,”
Mol. Cell. Biomech.
,
8
, pp.
215
232
.
42.
Su
,
T.
, and
Purohit
,
P. K.
,
2012
, “
Semiflexible Filament Networks Viewed as Fluctuating Beam-Frames
,”
Soft Matter
,
8
(
17
), pp.
4664
4674
.
43.
Magid
,
E.
,
Soldea
,
O.
, and
Rivlin
,
E.
,
2007
, “
A Comparison of Gaussian and Mean Curvature Estimation Methods on Triangular Meshes of Range Image Data
,”
Comput. Vision Image Understanding
,
107
(
3
), pp.
139
159
.
44.
Chen
,
L.
, and
Rong
,
Y.
,
2010
, “
Digital Topological Method for Computing Genus and the Betti Numbers
,”
Topol. Its Appl.
,
157
(
12
), pp.
1931
1936
.
45.
Lipowsky
,
R.
, and
Girardet
,
M.
,
1990
, “
Shape Fluctuations of Polymerized or Solidlike Membranes
,”
Phys. Rev. Lett.
,
65
(
23
), p.
2893
.
46.
Min
,
K.
, and
Aluru
,
N. R.
,
2011
, “
Mechanical Properties of Graphene Under Shear Deformation
,”
Appl. Phys. Lett.
,
98
(
1
), p.
013113
.
47.
Landau
,
L.
, and
Lifshitz
,
E.
,
1986
,
Theory of Elasticity
,
3rd ed.
,
Butterworth-Heinemann
,
Boston, MA
.
48.
Mounet
,
N.
, and
Marzari
,
N.
,
2005
, “
First-Principles Determination of the Structural, Vibrational and Thermodynamic Properties of Diamond, Graphite, and Derivatives
,”
Phys. Rev. B
,
71
(
20
), p.
205214
.