A Fluctuating Elastic Plate Model Applied to Graphene

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.

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.

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 w_i, 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 k_B 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 [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

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 $〈∬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⁻⁶) 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.

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 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 $ε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 [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_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].

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.

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

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

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.

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