Mechanics of thermally fluctuating membranes

Besides having unique electronic properties, graphene is claimed to be the strongest material in nature. In the press release of the Nobel committee it is claimed that a hammock made of a squared meter of one-atom thick graphene could sustain the wight of a 4 kg cat. More practically important are applications of graphene like scaffolds and sensors which are crucially dependent on the mechanical strength. Meter-sized graphene is even being considered for the lightsails in the starshot project to reach the star alpha centaury. The predicted strength of graphene is based on its very large Young modulus which is, per atomic layer, much larger than that of steel. This reasoning however would apply to conventional thin plates but does not take into account the peculiar properties of graphene as a thermally fluctuating crystalline membrane. It was shown recently both experimentally and theoretically that thermal fluctuations lead to a dramatic reduction of the Young modulus and increase of the bending rigidity for micron-sized graphene samples in comparison with atomic scale values. This makes the use of the standard F\"oppl-von Karman elasticity (FvK) theory for thin plates not directly applicable to graphene and other single atomic layer membranes. This fact is important because the current interpretation of experimental results is based on the FvK theory. In particular, we show that the FvK-derived Schwerin equation, routinely used to derive the Young modulus from indentation experiments has to be essentially modified for graphene at room temperature and for micron sized samples. Based on scaling analysis and atomistic simulation we investigate the mechanics of graphene under transverse load up to breaking. We determine the limits of applicability of the FvK theory and provide quantitative estimates for the different regimes.


II. INTRODUCTION
The deflection of a thin plate under transverse point or uniform load is normally well described by the FvK equations. These form a set of two coupled partial differential equations reading: where h is the displacement in the direction perpendicular to the plane, i.e. the deflection, φ is the potential for the in-plane stress tensor, κ is the bending rigidity, Y the two-dimensional (2D) Young modulus and P the transverse pressure. The behavior of h as a function of P and system size L can be obtained from a scaling analysis of the FvK equations as follows. Eq. 2 implies that φ/L 4 ∼ Y h 2 /L 4 or φ ∼ Y h 2 . Then the second term in eq. 1 scales as φh/L 4 ∼ Y h 3 /L 4 and dominates over the first term (∼ κh/L 4 ) in the regime of pressures yielding h 2 >> κ/Y . For graphene, with κ ≃ 1.1 eV and Y ≃ 19.6 eV/Å 2 at room temperature 16,17 , this condition implies h >> 0.23Å and is normally fulfilled for a system of mesoscopic size (or beyond), except for very low pressure. Hence, apart from this very low pressure regime, the deflection behaves as: h ≃ L 4 P gY 1/3 (3) where g is a dimensionless number depending on the shape of the 2D system and on the type of load, for instance uniform pressure or point load with a tip as in nano-indentation 2 . In the latter case, the pressure P in eq. 3 should be replaced by 4F/(πL 2 ) with F the force exerted by the tip. Then, eq. 3 turns into an equation which is equivalent to the so-called Schwerin equation, but without prestress. The Schwerin equation is commonly used to determine the elastic modulus Y from nano-indentation measurements 2,13 . From now on we will consider a circular drum of radius R = L/2, clamped at the edge, with h being the midpoint deflection. In that case, the value of g for point load has been derived from analytical solutions of the governing equations and is given by g = 16/(πg 3 ν ), withg ν ≃ 1.0491 − 0.1462ν − 0.1583ν 2 and ν the Poisson ratio 18 . With ν ≃ 0.26 for graphene, one findsg ν ≃ 1.0004 and g ≃ 5.087. For uniform load, a fully analytical solution is not available, but there are various approximate, semi-analytical solutions. The solutions reported in Refs. 8,18 yield g ν ≃ 0.7179 − 0.1706ν − 0.1495ν 2 ≃ 0.663 which is similar to the value obtained in Ref 19 yielding g ν ≃ 75(1 − ν 2 )/(8(23 + 18ν − 3ν 2 )) ≃ 0.686. It was noted 19 , however, that these expressions underestimate by 10 % the values obtained from the classical, more complex and accurate solution by Hencky 20 for the case that ν = 0.16. Here we will use g ν = 0.714 yielding g ≃ 43.9. With this larger value for g ν , compatible with Hencky's solution, the simulation data for h versus P yield the correct (known) value of Y which agrees with that obtained from a simulation with point load. We will comment on this later on.
At low enough pressure, where h 2 << κ/Y , we have a linear regime where h ≃ L 4 P/(f κ) according to the above scaling analysis with f another dimensionless, numerical factor. From the work in Ref. 8,18 , identifying the bending stiffness Y 3D d 3 /(12(1 − ν 2 )) for a thin plate with bulk Young modulus Y 3D and thickness d as κ for a membrane of atomic thickness, we can evaluate f ≃ 1024.0 and f = 256.0 for uniform and point load respectively.
Summing up the two terms from the scaling analysis, one finds that h(P, L) satisfies the equation: in agreement with the analysis in Refs. 18,21 . If we define the cross-over pressure P c1 as the pressure for which the linear term equals the non-linear term in eq. 4, which is the case when h 2 = f κ/(gY ), we find that the linear regime vanishes rapidly with system size as P c1 ≃ 2 (f κ) 3 /(gY )/L 4 . The above analysis is based on the assumption that the elastic moduli κ and Y are constant, i.e. independent of the system size. Recently it has been clearly confirmed 16 , however, that the elastic moduli of graphene are not material constants but scale as power-laws of the system size due to strong anharmonic coupling between in-plane modes and large out-of-plane modes, as predicted by membrane theory 22 . Besides, the moduli exhibit an anomalously strong dependences on strain 16 . Thus, generally speaking, for a 2D thermally fluctuating solid the above analysis is invalid and has to be adapted. Eventually this will lead to an anomalous deflection versus load relation h ∼ P α with α different from 1/3 (eq.3), as we will show explicitly below. Besides analytical results based on a revised scaling analysis for membranes, we also present the results from atomistic simulations for a graphene drum under uniform load, to validate our analytical findings.

III. RESULTS: SCALING THEORY
In order to account for the size and strain dependence of the elastic moduli, we extend our scaling analysis by replacing these moduli by their renormalized values κ R and Y R . The latter is given by 11,22 : while Y R /Y = 1 for L < L G and L σ < L G with η u ≃ 0.325 16 , where L G is the so-called Ginzburg length beyond which the power-law scaling is applicable. The length L σ is the size beyond which anharmonicity is suppressed due to tensile strain and is given by 23 : where ǫ is the average strain and where we used ). An equation similar to eq. 5 applies to κ R /κ, but with −η u replaced by η ≃ 1 − η u /2 ≃ 0.8375, implying that κ R increases with size while Y R decreases with size.
A theoretical estimate for L G , L theor 22 , yields L G ∼ 40Å at room temperature, but from simulations 16 , L G turned out to be about a factor 2 smaller. Therefore, in the further analysis we will use L G = c G L theor G , where c G ≃ 0.415 at 300 K is a correction factor resulting from analysis of the simulation data for Y R as a function of strain reported in Ref. 16 (see Supplementary Information S1).
With renormalized elastic constants, eq. 4 still holds, but with κ and Y replaced by κ R and Y R . Then, we can again determine the cross-over pressure P c1 imposing equality of the two terms on the left-hand size. For small load where L σ > L, applying the first line of eq. 5, one finds P c1 ∼ L −(6−η)/2 ∼ L −2.58 . For larger loads yielding L G < L σ < L, however, P c1 acquires a different size dependence: where The cross-over in the size dependence of P c1 should occur at a system size L c1 at another critical pressure, P c2 where L σ is equal to L c1 . Explicit expressions for P c2 as a function of L will be given below. For graphene, it turns out that L c1 would be smaller than L G , thus in a regime where renormalization does not apply. For L > L G > L c1 , P c2 < P c1 , implying that L σ < L at P c1 . Therefore, for graphene P c1 is always given by eq. 7, derived from the second line of eq. 5, and has no cross-over. Thus, with renormalized elastic moduli, the regime where the first term in eq. 4 is dominant still vanishes for L → ∞ but more slowly, namely as P c1 ∼ L −2.78 instead of L −4 .
A third critical pressure P c3 is defined as the pressure for which anharmonicity is completely suppressed, i.e. where L σ ≤ L G yielding Y R = Y . The important observation to make now is that for pressures P within P c2 < P < P c3 or equivalently L G < L σ < L, h as a function of P obeys a power-law different from that in eq. 3, due to the renormalization of the elastic moduli. Indeed, using eqs. 5 and 6, Y R depends on the strain ǫ as: In a similar way one can derive an expression for κ R (ǫ). Substitution of Eq. 8 with ǫ ≃ g ǫ h 2 /L 2 into Eq. 3 with Y replaced by Y R gives a self-consistency equation for h with solution: This equation replaces eq. 3 for the case of a 2D solid exhibiting renormalization of the elastic moduli according to membrane theory. One should notice that now the relation between h and P involves, apart from Y , also k B T /κ, which is natural as this quantity controls the strength of anharmonic coupling. Notice that eq. 3 is recovered from eq. 9 for µ = 0 (i.e. for η u = 0). The geometrical prefactor g ǫ for the strain in ǫ = g ǫ h 2 /L 2 depends on the shape of the deflected membrane. If we define ǫ = A/A 0 − 1 with A the surface area of the deflected membrane and A 0 = πL 2 /4 that of the flat drum, then for uniform load, with the shape of the drum approximately being that of a spherical cap, g ǫ = 2, whereas for nano-indentation, with an approximately cone shaped membrane, g ǫ = 1.
To find the critical pressure P c2 , we first need to solve L σ (h) = L for h. Substitution of this h into eq. 4, retaining both terms on the left-hand side with κ and Y replaced by κ R and Y R respectively, leads to: with (8 − 3η)/2 ≃ 2.74, (12 − 7η)/2 ≃ 3.07, g 21 ≃ 0.1063f (c 3η G f ν g ǫ ) −1/2 and g 22 ≃ 12.06g(c 7η−4 G f 3 ν g 3 ǫ ) −1/2 . In a similar way, for deriving P c3 we first solve L σ (h) = L G for h and then substitute this h into eq. 4. In this pressure regime the contribution from the term with κ is very small (for L not too small) and can be neglected. Then, we obtain: with g 3 ≃ 0.0146g(c 2 G f ν g ǫ ) −3/2 . If we just consider the pressure contribution ρmg due to the mass of the carbon atoms, with ρ the 2D density of graphene and g the gravitational acceleration, it can directly be calculated from eq. 11 that it requires a system size of about L = 305 km (!) to suppress anharmonicities by graphene's own weight.
The behavior for the various critical loads as a function of system size is depicted in Fig. 1 on two different length scales, corresponding to the scale used in our simulations and the typical scale in experiments respectively. For the latter case we used the parameters for point load and displayed critical forces F ci = πL 2 P ci /4 (i = 1, 2, 3) instead of pressures.

FIG. 2.
Snapshots from simulations of the indentation of a graphene drum with diameter L ≃ 315Å under uniform load (upper three graphs) and point load (lower three graphs). In the graphs with a view at an angle, one can see the thermal corrugation of the drum, responsible for the anharmonic effects. The bottom graphs show configurations just after breaking. For uniform load breaking occurs at the drum edge at P ≃ 60 kbar while for point load it occurs around the tip at F ≃ 160 nN for the given system size.

IV. RESULTS:ATOMISTIC SIMULATIONS
In order to validate the behavior derived above and in particular eqs. 9 we have performed atomistic, Monte Carlo simulations for a graphene drum with a diameter of L ≃ 315Å under uniform transverse pressures over a wide range between zero and 60 kbar. The LCBOPII model was used for the carbon interatomic interactions 24 and the pressure was modelled by assigning a weight M = P/(ρg) to each atom. Defining the z−direction perpendicular to the drum, a change δz of the z-coordinate of an atom contributes an amount M gδz to the energy change ∆E of the system entering the MC acceptance probability P acc = min[1, exp(−∆E/k B T )] for the configurational change. The simulations were performed at T = 300 K and the 2D density of the drum was adjusted to the equilibrium density at 300 K, so that no prestrain was present. For illustration, snapshots from these simulations are shown in Fig. 2, together with snapshots from a simulation under point load. In the latter case only atoms in a small circular, central region were assigned a weight M = F/(N c g), with F the total applied force and N c the number of atoms in the central circle. For our simulation, N c = 25. The simulation results for uniform pressure P versus h are given in Fig. 3. Before analyzing these data we should realize that for the derivation of eq. 9 we have tacitly neglected the linear term in eq. 4, which is justified for system sizes commonly used in experiments of the order of 1 µm or larger, yielding P c1 < 0.003 bar. The system size L ≃ 315 A used in our simulations, however, requires to include the linear regime to cover the pressure range smaller then P c1 ≃ 12 bar in this case. While for P < P c2 , P as a function of h should behave as P = Ah + Bh 3 , for P > P c2 the expected behavior is P = Ah (2−3η)/(2−η) + Bh 3+2µ ≃ Ah −0.441 + Bh 3.559 . In order to fit our simulation data for pressures below and above P c2 we used the following form which combines the usual FvK linear term with the renormalized expressions yielding correct asymptotic for h → 0 and h → ∞: It has two fitting parameters A and B of which the latter is related to the elastic moduli by B = (κ/k B T ) µg Y /L 4+2µ . The upper panel in Fig. 3 shows that for pressures up to ∼400 bar, the simulation data are in good agreement with eq. 12 as shown by the best fit (solid line), and clearly deviate from a best fit based on P = Ah + Bh 3 expected without renormalization of the elastic constants. Instead, beyond 400 bar up to about 6000 bar shown in the left bottom panel, the data points are best fitted by P = Ah + Bh 3 ≃ Bh 3 , which according to our analysis should apply for P > P c3 . This suggests that P c3 ≃ 400 bar, about a factor 2 smaller than the estimate from eq. 11 (see also Fig.  1) but nevertheless of the right order of magnitude. The numerical discrepancy here might be due to the fact that cross-over regimes are ignored in the theoretical derivations.
The best fit for pressures in the range [400, 40000] bar, beyond P c3 , based on P = Ah + Bh 3 yields B ≃ 0.128, implying a (bare) Young elastic modulus Y = L 4 B/g ≃ 299 N/m. This is the value after the mentioned adjustment of g ν such that it was equal to the Y obtained from a point load simulation at an applied force beyond F c3 . The found value is close to the known true bare modulus Y = 314 N/m at 300 K for LCBOPII. From the best fit at low pressures (< P c3 ) based on eq. 12 and with κ = 1.1 eV we obtain Y = 307 N/m. The small difference with the value from the first method may be due an uncertainty in the factor c g . In fact, the above two ways to determine Y (for known κ) can alternatively be used to determineg and by this c G (and L G ) by imposing equality of Y .

V. DISCUSSION
A way the extract the renormalization of Y from the simulation data is to make a fit based on eq. 12 and then use that Bh 3+2µ = Y R h 3 /(gL 4 ) to obtain the strain dependent Y R : valid for ǫ < ǫ c3 , where ǫ c3 ≃ 0.005 is the strain required to suppress anharmonicity, i.e. the strain beyond which the normal P ∼ h 3 is applicable. Notice that this approach does not require knowledge of κ norg. For point load, the same approach can be used, but with an additional factor π/(4L 2 ) multiplying the right-hand side of eq. 13.
The right bottom panel shows a deviation from the P ∼ h 3 behavior for deflections beyond 40Å. This deviation for large strain can be attributed to a normal softening of the elastic moduli due to stretching anharmonicity, as in 3D crystals. Such a softening for graphene under large strain is in agreement with previous observations 25 . The corresponding critical pressure, P c4 , should depend on the size as P c4 = g 4 /L, with g 4 ≃ 25.2 N/m derived from the value P c4 ≃ 8 kbar (8 10 8 N) for L = 3.15 10 −8 m. Indeed, assuming that this anharmonicity sets in at a fixed (size independent) critical strain, ǫ c4 = g ǫ h 2 c4 /L 2 , the size dependence of P c4 follows directly from eq. 4, neglecting the linear term, yielding P c4 ≃ (gY /L)(ǫ c4 /g ǫ ) 3/2 , with ǫ c4 ≃ 0.032. Although P c4 is a safe lower bound for breaking, normally breaking is only expected at significantly higher strains, where the in-plane moduli start to vanish due to the anharmonicity of the interaction potential. For LCBOPII, the bulk modulus vanishes at a strain value of ∼0.2, a value indeed close to the strain where breaking was actually observed in our simulation, at a pressure P br ≃ 50 kbar. This value of the breaking strain is similar to that found in a simulation study of graphene nanoribbons under uniaxial strain 26 . It should be noticed, however, that a typical atomistic simulation only covers a very small time interval, typically orders of magnitude smaller than a second, which makes the choice of a maximal, safe lower bound for breaking from simulations at a given temperature not obvious and somewhat arbitrary.
Staying on the safe side by choosing the breaking pressure as P br = 4P c4 = 4g 4 /L, yielding P br ≃ 32 kbar for the simulated system size with a corresponding strain of ∼0.13, a graphene drum of 1 m in diameter gives a breaking force of F c4 = πL 2 P c4 /4 ≃ 79.2 N, enough for an extremely heavy cat of about 8.0 kg to be safe, treating it as a uniform load. Treating the cat as a point load, however, and assuming that the breaking strain is equal to that for uniform load, we have to correct g 4 by a factor ∼ 0.328 due to the different values for g and g ǫ , implying that the cat should not be heavier than 2.65 kg, i.e. a young cat, to be safe. In reality, a cat on a drum of this size is something between point and uniform load, so that probably any cat should be safe on it. It is interesting to notice, and somewhat counterintuitive, that while a drum of 1 m cannot bear a person of 100 kg, a drum of 40 m could, due to the fact that F c4 grows linearly with L. A graphene drum would only break by its own weight for a size L = 4g 4 /(ρmg) ≃ 13520 km ! While the relations derived above are appropriate for the analysis of our simulations where prestress can be controlled and taken to be zero, it should be noticed that in nano-indentation experiments almost unavoidably some prestress σ 0 is present, created during preparation of the drum. The implications for the load versus deflection expression and the various critical loads for the case of tensile prestress, including renormalization of the elastic moduli, are given in the Supplementary Information S3. Tensile prestress gives rise to a contribution to the force which is linear in h, namely πσ 0 h. As this is an order L 2 larger then the linear contribution κh/L 2 arising from the FvK equations, its contribution increases with system size as L 2η−1 , and can be significant as compared to the cubic term Y h 3 /L 2 , even for µm sized drums. Therefore, for the sake of accuracy in measuring the elastic modulus and the effect of its renormalization, one should keep the prestress as small as possible, so that the term cubic in h, from which the elastic modulus is determined, is the dominant term. Moreover, while F c1 increases with system size for tensile prestress, F c2 and F c3 decrease with σ 0 . To be able to measure the renormalization of the elastic modulus, however, F c3 should not be too small, leaving a sufficiently large force domain (< F c3 ) for observing renormalization.

VI. CONCLUSIONS
The revised FvK theory for thermally excited membranes like graphene that we have presented here is important for any technological application of 2D materials involving their mechanical properties. For graphene, we have discussed in a quantitative way the behaviour under uniform and point load up to breaking.