A Fluctuating Elastic Plate Model Applied to Graphene

Xiaojun Liang; Prashant K. Purohit

doi:10.1115/1.4033681

Journal of Applied Mechanics | Volume 83 | Issue 8 | research-article

< Previous Article Next Article >

Research Papers

A Fluctuating Elastic Plate Model Applied to Graphene OPEN ACCESS

Xiaojun Liang and Prashant K. Purohit

[+] Author and Article Information

Xiaojun Liang

Department of Mechanical Engineering
and Applied Mechanics,
University of Pennsylvania,
Philadelphia, PA 19104

Prashant K. Purohit

Department of Mechanical Engineering
and Applied Mechanics,
University of Pennsylvania,
Philadelphia, PA 19104
e-mail: purohit@seas.upenn.edu

¹Corresponding author.

Contributed by the Applied Mechanics Division of ASME for publication in the JOURNAL OF APPLIED MECHANICS. Manuscript received March 24, 2016; final manuscript received May 18, 2016; published online June 9, 2016. Editor: Yonggang Huang.

J. Appl. Mech 83(8), 081008 (Jun 09, 2016) (8 pages) Paper No: JAM-16-1155; doi: 10.1115/1.4033681 History: Received March 24, 2016; Revised May 18, 2016

Article

References

Figures

Tables

Citing Articles

Abstract

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.

FIGURES IN THIS ARTICLE

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.

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: Display Formula

\begin{matrix} E_{von} = E_{s} + E_{b} \\ = \iint d x d y \frac{Y}{2} [{(ε_{x} + ε_{y})}^{2} - 2 (1 - ν) (ε_{x} ε_{y} - ε_{x y}^{2})] \\ + \iint d x d y [\frac{K_{b}}{2} {(w_{, x x} + w_{, y y})}^{2} + K_{G} (w_{, x x} w_{, y y} - w_{, x y}^{2})] \end{matrix}

where Y and ν are the in-plane Young's modulus and Poisson's ratio of the material, respectively, K_b is a mean curvature bending modulus, and K_G is the Gauss curvature bending modulus. E_s is the energy of in-plane strains, and E_b is the energy of out-of-plane bending. The variables $ε_{x}, ε_{y}, ε_{x y}$ , and w in the integral are, respectively, the in-plane strains and out-plane displacement of the neutral plane. The subindices $_{, x x},_{, x y},_{, y y}$ represent second derivatives, respectively, with respect to the reference coordinates x and y. We recall that the strains ${ε_{x}, ε_{y}, ε_{x y}}$ are not independent since Display Formula

\begin{matrix} ε_{x} = u_{, x} + \frac{w_{, x}^{2}}{2} \\ ε_{y} = v_{, y} + \frac{w_{, y}^{2}}{2} \\ ε_{x y} = \frac{u_{, y} + v_{, x}}{2} + \frac{w_{, x} w_{, y}}{2} \end{matrix}

These strain–displacement relations are connected through the compatibility equation Display Formula

ε_{x, y y} + ε_{y, x x} - 2 ε_{x y, x y} = w_{, x y}^{2} - w_{, x x} w_{, y y}

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 Display Formula

\begin{matrix} E_{f} = - \iint F (u_{, x} + v_{, y}) d x d y \\ = - \iint F (ε_{x} - \frac{w_{, x}^{2}}{2} + ε_{y} - \frac{w_{, y}^{2}}{2}) d x d y \end{matrix}

Equations (1) and (4) together give the energy of a configuration of the plate which is determined by $ε_{x} (x, y), ε_{y} (x, y), ε_{x y} (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 $E_{s} \sim Y A ε^{2} / 2$ , where ε represents in-plane strains, and $E_{b} \sim K_{b} A κ^{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} 〉 = k_{B} T / Y A$ and $〈 κ^{2} 〉 = k_{B} T / K_{b} A$ . 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} 〉 h^{2} = K_{b} / Y h^{2}$ . For a graphene sheet, Y ≈ 1 TPa × 0.3 nm, K_b ≈ 10⁻¹⁹ N·m, and h ≈ 0.3 nm [28], so that $K_{b} / Y h^{2} \approx 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 Display Formula

\begin{matrix} E = E_{b} + E_{f} + C = \iint d x d y [\frac{K_{b}}{2} {(w_{, x x} + w_{, y y})}^{2} + K_{G} (w_{, x x} w_{, y y} - w_{, x y}^{2}) + F (\frac{w_{, x}^{2}}{2} + \frac{w_{, y}^{2}}{2})] + C \end{matrix}

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 $4 N^{2} / \sqrt{3}$ equilateral triangles, as shown in Fig. 1(a), with $P \approx 2 N (N + 1) / \sqrt{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 = w M w^{T}$ , where the vector $w = [w_{1}, w_{2}, \dots, w_{P}]$ contains all the node displacements. The matrix M is the stiffness matrix, and it is a function of $K_{b}, K_{G}, F, L, and l$ . The probability of finding the membrane in a given configuration of energy E is proportional to $exp (- E / k_{B} T)$ , where k_B is the Boltzmann constant and T is the absolute temperature. The partition function Z is obtained by integration of $exp (- E / k_{B} T)$ 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

Display Formula

Z = \sqrt{\frac{{(2 π k_{B} T)}^{P}}{\det M}}

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 Display Formula

f_{c} = ε_{x, y y} + ε_{y, x x} - 2 ε_{x y, x y} - w_{, x y}^{2} + w_{, x x} w_{, y y}

The energy cost can be computed by multiplying 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 Display Formula

\bar{E} = E + \iint λ f_{c} d A = E + \sum_{i} \iint λ_{i} f_{c}^{i} d A_{i}

where λ(x, y) is a scalar with energy units, and A_i 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 A_i such that it could come out of each integral. Then, we express $\iint f_{c}^{i} d A_{i}$ as a quadratic expression in terms of the displacement vector w just as we have done for the plate energy EDisplay Formula

\bar{E} = w M w^{T} + \sum_{i} λ_{i} w_{i} M_{i}^{c} {w_{i}}^{T}

Here, $w_{i} = [w, w_{1}, w_{2}, w_{3}, w_{4}, w_{5}, w_{6}]$ is a subset of w which consists of nodal displacement of node i and six adjacent nodes as shown in Fig. 1(b). $M_{i}^{c}$ is a stiffness matrix of compatibility at node i, such that $\iint f_{c}^{i} d A_{i} = w_{i} M_{i}^{c} {w_{i}}^{T}$ . $M_{i}^{c}$ 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 $w_{i} M_{i}^{c} {w_{i}}^{T}$ . This is done by taking the partial derivative of the Gibb's free energy with respect to λ_i as in Refs. [41] and [42] Display Formula

〈 w_{i} M_{i}^{c} {w_{i}}^{T} 〉 = \frac{\partial G}{\partial λ_{i}}

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 Display Formula

\iint C_{G} d A = 2 π - \sum γ_{i}

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, w₁, and w₂ 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₁, and w₂ as Display Formula

\begin{matrix} a = \sqrt{{(w - w_{1})}^{2} + l^{2}} \\ b = \sqrt{{(w - w_{2})}^{2} + l^{2}} \\ c = \sqrt{{(w_{1} - w_{2})}^{2} + l^{2}} \end{matrix}

So, the deformed angle is given by Display Formula

γ = \arccos \frac{a^{2} + b^{2} - c^{2}}{2 a b}

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 Display Formula

γ - \frac{π}{3} = \frac{2 {(w_{1} - w_{2})}^{2} - {(w - w_{1})}^{2} - {(w - w_{2})}^{2}}{2 \sqrt{3} l^{2}}

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 Display Formula

\begin{matrix} \iint C_{G} d A = \sum_{i = 1}^{6} (\frac{π}{3} - γ_{i}) \\ = \frac{\sum_{i = 1}^{6} {(w - w_{i})}^{2} - \sum_{i = 1}^{6} {(w_{i} - w_{i + 1})}^{2}}{\sqrt{3} l^{2}} \end{matrix}

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 Display Formula

\begin{matrix} I = \iint (ε_{x, y y} + ε_{y, x x} - 2 ε_{x y, x y}) d A \\ = \oint (ε_{x y, x} - ε_{x, y}) d x + \oint (ε_{y, x} - ε_{x y, y}) d y \end{matrix}

The terms $ε_{x y, x}$ and $ε_{x y, y}$ in the integrands will contribute zero because our integration contour is closed. Hence, we are left with Display Formula

I = \oint ε_{y, x} d y - \oint ε_{x, y} d x

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 Display Formula

I = \oint [\frac{1}{2} \frac{\partial}{\partial x} {(\frac{\partial w}{\partial y})}^{2} d y - \frac{1}{2} \frac{\partial}{\partial y} {(\frac{\partial w}{\partial x})}^{2} d x]

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 f_c, could be evaluated up to quadratic order in the displacements $w_{i}$ at the nodes, governed by matrix $M_{i}^{c}$ 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 Display Formula

\bar{E} = w M w^{T} + λ \sum_{i} w_{i} M_{i}^{c} {w_{i}}^{T}

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 Display Formula

\sum_{i} w_{i} M_{i}^{c} {w_{i}}^{T} = w M_{c} w^{T}

where the matrix $M_{c}$ has nonzero entries only for the boundary nodes. Furthermore, we find that $M_{c}$ 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 $〈 \iint f_{c}^{i} d A_{i} 〉$ 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⁻⁶) 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_c 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 f_c over all nodes approaches zero as λ becomes sufficiently large (see the dashed line associated with right axis). We see that λ = 10⁹ 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 $〈 w_{i} w_{j} 〉$ [25,40–42], which can be calculated from the inverse of the stiffness matrix M as Display Formula

〈 w_{i} w_{j} 〉 = \frac{k_{B} T}{{[M^{- 1}]}_{i j}}

From this, we can get the variance of the out-plane deflection $〈 w_{i}^{2} 〉$ 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_G in each plot for comparison. Recall that since a fluid membrane has zero in-plane shear modulus, we can choose $ε_{x y} (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 $\bar{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 [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].

In the analysis given by Boal et al. [5,45], an ideal membrane with nonzero in-plane shear resistance μ, but zero bending resistance K_b, has been shown to have an average thermal fluctuation $\bar{w} \propto L^{0.5}$ given by Display Formula

\bar{w} = 0.62 L^{0.5} {(\frac{k_{B} T}{Y})}^{0.25}

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_b ≠ 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 K_b and Gaussian curvature modulus K_G to be exactly zero since this results in a singular stiffness matrix M. So, we choose K_b and K_G to be smaller than the thermal energy scale 1 k_BT. 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].

Study of Ripples Under Shear Loading.

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 Display Formula

\begin{matrix} E_{f} = - \iint S (u_{, x} - v_{, y}) d x d y \\ = - \iint S (ε_{x} - \frac{w_{, x}^{2}}{2} - ε_{y} + \frac{w_{, y}^{2}}{2}) d x d y \end{matrix}

This results in a sign change in Eq. (5) as Display Formula

\begin{matrix} E = \iint d x d y [\frac{K_{b}}{2} {(w_{, x x} + w_{, y y})}^{2} \\ + K_{G} (w_{, x x} w_{, y y} - w_{, x y}^{2}) + S (\frac{w_{, x}^{2}}{2} - \frac{w_{, y}^{2}}{2})] \end{matrix}

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]

Display Formula

w (x, y) = \sum_{(m, n)} q_{m n} sin (\frac{n π x}{L_{1}}) sin (\frac{m π y}{L_{2}})

where q_mn is the amplitude of the normal mode corresponding to m, n. Note that we assume our plate to be rectangular with sides L₁ in x direction and L₂ in y direction. When we plug the above expression into Eq. (24), we get Display Formula

E = \sum_{(m, n)} q_{m n}^{2} [\frac{K b π^{4}}{8 L_{1} L_{2}} {(n^{2} + m^{2})}^{2} + \frac{S π^{2}}{8} (n^{2} - m^{2})]

From the equipartition theorem of statistical mechanics [47], each Fourier mode (m, n) should possess energy k_BT/2. Therefore, we can compute the amplitude of each mode as Display Formula

〈 q_{m n}^{2} 〉 = \frac{4 k_{B} T L_{1} L_{2}}{K_{b} π^{4} {(n^{2} + m^{2})}^{2} + S L_{1} L_{2} π^{2} (n^{2} - m^{2})}

Therefore, the mode (m, n) with maximum amplitude should be given by Display Formula

n = 1; m = [\frac{\sqrt{1 + \frac{2 S L_{1} L_{2}}{K_{b} π^{2}}} - 1}{2}]

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⁻⁶ K⁻¹. 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 = - k_{B} T \ln Z$ . 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 (\infty, 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 Display Formula

\frac{Δ A}{A} = - \frac{k_{B} T}{2 A} \frac{\partial}{\partial F} \ln (\det M)

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 K_b 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: Display Formula

K_{b} = 131.36 + 200 \cdot tan h \frac{T}{1500}

where the units of K_b 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 Display Formula

α = - \frac{\partial^{2}}{\partial T \partial F} [\frac{k_{B} T}{2 A} \ln det M]

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⁻¹, 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∂²G/∂T² at given tension F or shear loading S through

Display Formula

\frac{\partial α}{\partial T} = - \frac{1}{A} \frac{\partial^{3} G}{\partial F \partial T^{2}} = \frac{1}{A T} \frac{\partial C_{F}}{\partial F}

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.

Acknowledgements

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

References

Helfrich, W. , 1973, “ Elastic Properties of Lipid Bilayers: Theory and Possible Experiments,” Z. Naturforsch. C, 28(11), pp. 693–703. [PubMed]

Weeks, J. D. , 1977, “ Structure and Thermodynamics of the Liquid–Vapor Interface,” J. Chem. Phys., 67(7), pp. 3106–3121. [CrossRef]

Evans, E. , and Rawicz, W. , 1990, “ Entropy-Driven Tension and Bending Elasticity in Condensed-Fluid Membranes,” Phys. Rev. Lett., 64(17), p. 2094. [CrossRef] [PubMed]

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. [CrossRef] [PubMed]

Boal, D. , 2002, Mechanics of the Cell, Cambridge University Press, Cambridge, UK.

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.

Kosmrlj, A. , and Nelson, D. R. , 2014, “ Thermal Excitations of Warped Membranes,” Phys. Rev. E, 89(2), p. 022126. [CrossRef]

Paulose, J. , Vliegenthart, G. A. , Gompper, G. , and Nelson, D. R. , 2012, “ Fluctuating Shells Under Pressure,” Proc. Natl. Acad. Sci., 109(48), pp. 19551–19556. [CrossRef]

Helfrich, W. , 1975, “ Out-of-Plane Fluctuations of Lipid Bilayers,” Z. Naturforsch. Sect. C: Biosci., 30(6), p. 841.

Milner, S. T. , and Safran, S. A. , 1987, “ Dynamical Fluctuations of Droplet Microemulsions and Vesicles,” Phys. Rev. A, 36(9), p. 4371. [CrossRef]

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. [CrossRef]

Auth, T. , and Gompper, G. , 2013, “ Fluctuation Pressure of Biomembranes in Planar Confinement,” Phys. Rev. E, 88(1), p. 010701. [CrossRef]

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. [CrossRef]

Zhang, T. , Li, X. , and Gao, H. , 2014, “ Defects Controlled Wrinkling and Topological Design in Graphene,” J. Mech. Phys. Solids, 67, pp. 2–13. [CrossRef]

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. [CrossRef] [PubMed]

Ramakrishnan, N. , Kumar, P. S. , and Radhakrishnan, R. , 2014, “ Mesoscale Computational Studies of Membrane Bilayer Remodeling by Curvature-Inducing Proteins,” Phys. Rep., 543(1), pp. 1–60. [CrossRef] [PubMed]

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. [CrossRef]

Hanlumyuang, Y. , Liu, L. , and Sharma, P. , 2014, “ Revisiting the Entropic Force Between Fluctuating Biological Membranes,” J. Mech. Phys. Solids, 63, pp. 179–186. [CrossRef]

Audoly, B. , and Pomeau, Y. , 2010, Elasticity and Geometry, Oxford University Press, Oxford, UK.

Nelson, D. R. , and Peliti, L. , 1987, “ Fluctuations in Membranes With Crystalline and Hexatic Order,” J. Phys., 48(7), pp. 1085–1092. [CrossRef]

Fasolino, A. , Los, J. H. , and Katsnelson, M. I. , 2007, “ Intrinsic Ripples in Graphene,” Nat. Mater., 6(11), pp. 858–861. [CrossRef] [PubMed]

Le Doussal, P. , and Radzihovsky, L. , 1992, “ Self-Consistent Theory of Polymerized Membranes,” Phys. Rev. Lett., 69(8), p. 1209. [CrossRef] [PubMed]

Abraham, F. F. , and Nelson, D. R. , 1990, “ Diffraction From Polymerized Membranes,” Science, 249(4967), pp. 393–397. [CrossRef] [PubMed]

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. [CrossRef]

Su, T. , and Purohit, P. K. , 2010, “ Thermomechanics of a Heterogeneous Fluctuating Chain,” J. Mech. Phys. Solids, 58(2), pp. 164–186. [CrossRef]

Agrawal, N. J. , and Radhakrishnan, R. , 2009, “ Calculation of Free Energies in Fluid Membranes Subject to Heterogeneous Curvature Fields,” Phys. Rev. E, 80(1), p. 011925. [CrossRef]

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. [CrossRef] [PubMed]

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. [CrossRef] [PubMed]

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. [CrossRef]

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. [CrossRef]

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. [CrossRef] [PubMed]

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. [CrossRef]

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. [CrossRef] [PubMed]

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.

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. [CrossRef]

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. [CrossRef] [PubMed]

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.

Fraternali, F. , and Marcelli, G. , 2012, “ A Multiscale Approach to the Elastic Moduli of Biomembrane Networks,” Biomech. Model. Mechanobiol., 11(7), pp. 1097–1108. [CrossRef] [PubMed]

Flory, P. J. , and Volkenstein, M. , 1969, “ Statistical Mechanics of Chain Molecules,” Biopolymers, 8(5), pp. 699–700. [CrossRef]

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. [CrossRef] [PubMed]

Su, T. , and Purohit, P. K. , 2011, “ Fluctuating Elastic Filaments Under Distributed Loads,” Mol. Cell. Biomech., 8, pp. 215–232. [PubMed]

Su, T. , and Purohit, P. K. , 2012, “ Semiflexible Filament Networks Viewed as Fluctuating Beam-Frames,” Soft Matter, 8(17), pp. 4664–4674. [CrossRef]

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. [CrossRef]

Chen, L. , and Rong, Y. , 2010, “ Digital Topological Method for Computing Genus and the Betti Numbers,” Topol. Its Appl., 157(12), pp. 1931–1936. [CrossRef]

Lipowsky, R. , and Girardet, M. , 1990, “ Shape Fluctuations of Polymerized or Solidlike Membranes,” Phys. Rev. Lett., 65(23), p. 2893. [CrossRef] [PubMed]

Min, K. , and Aluru, N. R. , 2011, “ Mechanical Properties of Graphene Under Shear Deformation,” Appl. Phys. Lett., 98(1), p. 013113. [CrossRef]

Landau, L. , and Lifshitz, E. , 1986, Theory of Elasticity, 3rd ed., Butterworth-Heinemann, Boston, MA.

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. [CrossRef]

View article in PDF format.

View PDF