Mechanics of thermally fluctuating membranes

TL;DR

This paper investigates how thermal fluctuations affect the mechanical properties of graphene membranes, showing that standard thin plate theories need modification for accurate modeling at room temperature and micron scales.

Contribution

The study reveals the limitations of Föppl-von Kármán theory for graphene, providing modified equations and scaling analysis to better understand its mechanics under load.

Findings

Thermal fluctuations significantly reduce Young's modulus in micron-sized graphene.

Standard FvK theory does not accurately describe graphene mechanics at room temperature.

Modified models improve the interpretation of indentation experiments.

Abstract

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.

Equations

$$\kappa {\Delta }^{2}h-\left(\frac{{\partial }^2\varphi }{\partial y^2}\frac{{\partial }^{2}h}{\partial x^2}+\frac{{\partial }^2\varphi }{\partial x^2}\frac{{\partial }^{2}h}{\partial y^2}-\frac{{\partial }^2\varphi }{\partial x\partial y}\frac{{\partial }^{2}h}{\partial x\partial y}\right)=P$$

$${\Delta }^2\varphi -Y\left(\frac{{\partial }^{2}h}{\partial x^2}\frac{{\partial }^{2}h}{\partial y^2}-{\left(\frac{{\partial }^{2}h}{\partial x\partial y}\right)}^2\right)=0$$

$$h\simeq {\left(\frac{L^{4}P}{gY}\right)}^{1/3}$$

$$f\kappa h+gY{h}^3=L^{4}P$$

$$\frac{Y_{R}}{Y}\sim\left\{\begin{array}[]{lc}\vspace*{0.1cm}\displaystyle\left(\frac{L}{L_{G}}\right)^{-\eta_{u}}&~{}~{}L_{G}<L<L_{\sigma}\\ \vspace*{0.1cm}\displaystyle\left(\frac{L_{\sigma}}{L_{G}}\right)^{-\eta_{u}}&~{}~{}L_{G}<L_{\sigma}\leq L\end{array}\right.$$

$$L_{\sigma }={\left(\frac{(2\pi {)}^2\kappa }{2BϵL_G^{\eta }}\right)}^{\frac{1}{2-\eta }}={\left(\frac{(2\pi {)}^2\kappa }{f_{\nu }YϵL_G^{\eta }}\right)}^{\frac{1}{2-\eta }}$$

$$P_{c1}=\frac{{\tilde{g}}_1{\kappa }^{\frac{8-8\eta }{8-4\eta }}{Y}^{\frac{3\eta -2}{8-4\eta }}(k_{B}T{)}^{\frac{2+\eta }{8-4\eta }}}{L^{\frac{14-9\eta }{4-2\eta }}}$$

$$Y_R(ϵ)\simeq Y{\left(\frac{(2\pi {)}^2\kappa }{f_{\nu }YϵL_G^2}\right)}^{-\mu }\simeq Y{\left(\frac{16\pi \kappa c_G^2{f}_{\nu }ϵ}{3k_{B}T}\right)}^{\mu }$$

$$h\simeq {\left(\frac{k_{B}T}{\kappa }\right)}^{\frac{\mu }{3+2\mu }}{\left(\frac{L^{4+2\mu }P}{\tilde{g}Y}\right)}^{\frac{1}{3+2\mu }}$$

$$P_{c2}=\frac{g_{21}(k_{B}T{)}^{\frac{3\eta }{4}}{\kappa }^{\frac{3-3\eta }{4}}{Y}^{\frac{3\eta -2}{4}}}{L^{\frac{8-3\eta }{2}}}+\frac{g_{22}(k_{B}T{)}^{\frac{7\eta -4}{4}}{\kappa }^{\frac{7-7\eta }{2}}}{Y^{\frac{6-7\eta }{4}}{L}^{\frac{12-7\eta }{2}}}$$

$$P_{c3}=\frac{g_{3}Y}{L}{\left(\frac{k_{B}T}{\kappa }\right)}^{3/2}$$

$$P=Ah+B{h}^{3+2\mu }$$

$$Y_R(ϵ)=\frac{L^{4+2\mu }B}{g}{\left(\frac{h^2}{L^2}\right)}^{\mu }=\frac{L^{4+2\mu }B}{g{g}_{ϵ}^{\mu }}{ϵ}^{\mu }$$

Full text

Mechanics of thermally fluctuating membranes.

J. H. Los, A. Fasolino∗, and M. I. Katsnelson

Radboud University, Institute for Molecules and Materials, Heyendaalseweg 135, 6525AJ Nijmegen, The Netherlands

I Abstract

Besides having unique electronic properties, graphene is claimed to be the strongest material in naturebooth2008 ; lee2008 . In the press release of the Nobel committeenobel it is claimed that a hammock made of a squared meter of one-atom thick graphene could sustain the weight of a 4 kg cat. More practically important are many applications of graphene like scaffoldsscaffold and sensorssensor which are crucially dependent on the mechanical strength. Meter-sized graphene is even being considered as material for the lightsails in the starshot project to reach the star alpha centaurystarshot . The predicted exceptional 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 platesLandauelasticity ; Timoscenko but does not take into account the peculiar properties of graphene as a thermally fluctuating crystalline membranebookMIsha ; accountchemreas ; nelsonribbons ; mirlin-arxiv-Hooks . It was shown recently both experimentallyspanishnatphys ; mcewuenkirigami ; bolotinNat and theoreticallyLos1 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öppl-von Karman elasticity (FvK) theory for thin platesLandauelasticity ; Timoscenko 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 experimentsspanishnatphys 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:

[TABLE]

[TABLE]

where $h$ is the displacement in the direction perpendicular to the plane, i.e. the deflection, $\phi$ is the potential for the in-plane stress tensor, $\kappa$ 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 $\phi/L^{4}\sim Yh^{2}/L^{4}$ or $\phi\sim Yh^{2}$. Then the second term in eq. 1 scales as $\phi h/L^{4}\sim Yh^{3}/L^{4}$ and dominates over the first term ($\sim\kappa h/L^{4}$) in the regime of pressures yielding $h^{2}>>\kappa/Y$. For graphene, with $\kappa\simeq 1.1$ eV and $Y\simeq 19.6$ eV/ Å2 at room temperatureFasolino1 ; Los1 , 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:

[TABLE]

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 lee2008 . In the latter case, the pressure $P$ in eq. 3 should be replaced by $4F/(\pi 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 measurementslee2008 ; spanishnatphys . 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/(\pi\tilde{g}_{\nu}^{3})$, with $\tilde{g}_{\nu}\simeq 1.0491-0.1462\nu-0.1583\nu^{2}$ and $\nu$ the Poisson ratioKomaragiri . With $\nu\simeq 0.26$ for graphene, one finds $\tilde{g}_{\nu}\simeq 1.0004$ and $g\simeq 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. Timoscenko ; Komaragiri yield $g_{\nu}\simeq 0.7179-0.1706\nu-0.1495\nu^{2}\simeq 0.663$ which is similar to the value obtained in Ref Yue yielding $g_{\nu}\simeq 75(1-\nu^{2})/(8(23+18\nu-3\nu^{2}))\simeq 0.686$. It was noted Yue , however, that these expressions underestimate by 10 % the values obtained from the classical, more complex and accurate solution by Hencky Hencky for the case that $\nu=0.16$. Here we will use $g_{\nu}=0.714$ yielding $g\simeq 43.9$. With this larger value for $g_{\nu}$, 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}<<\kappa/Y$, we have a linear regime where $h\simeq L^{4}P/(f\kappa)$ according to the above scaling analysis with $f$ another dimensionless, numerical factor. From the work in Ref.Timoscenko ; Komaragiri , identifying the bending stiffness $Y_{3D}d^{3}/(12(1-\nu^{2}))$ for a thin plate with bulk Young modulus $Y_{3D}$ and thickness $d$ as $\kappa$ for a membrane of atomic thickness, we can evaluate $f\simeq 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:

[TABLE]

in agreement with the analysis in Refs.Komaragiri ; Wang . 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\kappa/(gY)$, we find that the linear regime vanishes rapidly with system size as $P_{c1}\simeq 2\sqrt{(f\kappa)^{3}/(gY)}/L^{4}$.

The above analysis is based on the assumption that the elastic moduli $\kappa$ and $Y$ are constant, i.e. independent of the system size. Recently it has been clearly confirmedLos1 , 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 Nelson1 . Besides, the moduli exhibit an anomalously strong dependences on strainLos1 . 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\sim P^{\alpha}$ with $\alpha$ 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 $\kappa_{R}$ and $Y_{R}$. The latter is given byNelson1 ; nelsonribbons :

[TABLE]

while $Y_{R}/Y=1$ for $L<L_{G}$ and $L_{\sigma}<L_{G}$ with $\eta_{u}\simeq 0.325$Los1 , where $L_{G}$ is the so-called Ginzburg length beyond which the power-law scaling is applicable. The length $L_{\sigma}$ is the size beyond which anharmonicity is suppressed due to tensile strain and is given byRoldan :

[TABLE]

where $\epsilon$ is the average strain and where we used the relation $2B=Y/(1-\nu)=f_{\nu}Y$ with $f_{\nu}\equiv 1/(1-\nu)\simeq 1.35$ (see also Supplementary Information S2). An equation similar to eq. 5 applies to $\kappa_{R}/\kappa$, but with $-\eta_{u}$ replaced by $\eta\simeq 1-\eta_{u}/2\simeq 0.8375$, implying that $\kappa_{R}$ increases with size while $Y_{R}$ decreases with size.

A theoretical estimate for $L_{G}$, $L_{G}^{theor}=2\pi\sqrt{16\pi\kappa^{2}/(3Yk_{B}T)}$Nelson1 , yields $L_{G}\sim$ 40 Å  at room temperature, but from simulationsLos1 , $L_{G}$ turned out to be about a factor 2 smaller. Therefore, in the further analysis we will use $L_{G}=c_{G}L_{G}^{theor}$, where $c_{G}\simeq 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.Los1 (see Supplementary Information S1).

With renormalized elastic constants, eq. 4 still holds, but with $\kappa$ and $Y$ replaced by $\kappa_{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_{\sigma}>L$, applying the first line of eq. 5, one finds $P_{c1}\sim L^{-(6-\eta)/2}\sim L^{-2.58}$. For larger loads yielding $L_{G}<L_{\sigma}<L$, however, $P_{c1}$ acquires a different size dependence:

[TABLE]

where $(14-9\eta)/(4-2\eta)\simeq 2.78$ and $\tilde{g}_{1}\simeq 0.55f^{\nu_{1}}g^{\nu_{2}}$ $(c_{G}^{2}g_{\sigma}g_{\epsilon})^{-\nu_{3}}$, with $\nu_{1}=(10-7\eta)/(8-4\eta)$, $\nu_{2}=(3\eta-2)/(8-4\eta)$ and $\nu_{3}=(2+\eta)/(8-4\eta)$. 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_{\sigma}$ 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_{\sigma}<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\rightarrow\infty$ but more slowly, namely as $P_{c1}\sim 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_{\sigma}\leq 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_{\sigma}<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 $\epsilon$ as:

[TABLE]

with $\mu\equiv(2-2\eta)/(2-\eta)\simeq 0.2797$. In a similar way one can derive an expression for $\kappa_{R}(\epsilon)$. Substitution of Eq. 8 with $\epsilon\simeq g_{\epsilon}h^{2}/L^{2}$ into Eq. 3 with $Y$ replaced by $Y_{R}$ gives a self-consistency equation for $h$ with solution:

[TABLE]

where $(3+2\mu)\simeq 3.56$ and $\tilde{g}\simeq(16\pi c_{G}^{2}f_{\nu}g_{\epsilon}/3)^{\mu}g$. 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/\kappa$, which is natural as this quantity controls the strength of anharmonic coupling. Notice that eq. 3 is recovered from eq. 9 for $\mu=0$ (i.e. for $\eta_{u}=0$).

The geometrical prefactor $g_{\epsilon}$ for the strain in $\epsilon=g_{\epsilon}h^{2}/L^{2}$ depends on the shape of the deflected membrane. If we define $\epsilon=\sqrt{A/A_{0}}-1$ with $A$ the surface area of the deflected membrane and $A_{0}=\pi 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_{\epsilon}=2$, whereas for nano-indentation, with an approximately cone shaped membrane, $g_{\epsilon}=1$.

To find the critical pressure $P_{c2}$, we first need to solve $L_{\sigma}(h)=L$ for $h$. Substitution of this $h$ into eq. 4, retaining both terms on the left-hand side with $\kappa$ and $Y$ replaced by $\kappa_{R}$ and $Y_{R}$ respectively, leads to:

[TABLE]

with $(8-3\eta)/2\simeq 2.74$, $(12-7\eta)/2\simeq 3.07$, $g_{21}\simeq 0.1063f(c_{G}^{3\eta}f_{\nu}g_{\epsilon})^{-1/2}$ and $g_{22}\simeq 12.06g(c_{G}^{7\eta-4}f_{\nu}^{3}g_{\epsilon}^{3})^{-1/2}$.

In a similar way, for deriving $P_{c3}$ we first solve $L_{\sigma}(h)=L_{G}$ for $h$ and then substitute this $h$ into eq. 4. In this pressure regime the contribution from the term with $\kappa$ is very small (for $L$ not too small) and can be neglected. Then, we obtain:

[TABLE]

with $g_{3}\simeq 0.0146g(c_{G}^{2}f_{\nu}g_{\epsilon})^{-3/2}$. If we just consider the pressure contribution $\rho mg$ due to the mass of the carbon atoms, with $\rho$ 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}=\pi L^{2}P_{ci}/4~{}(i=1,2,3)$ instead of pressures.

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\simeq 315$ Å  under uniform transverse pressures over a wide range between zero and 60 kbar. The LCBOPII model was used for the carbon interatomic interactionsLos2 and the pressure was modelled by assigning a weight $M=P/(\rho g)$ to each atom. Defining the $z-$direction perpendicular to the drum, a change $\delta z$ of the $z$-coordinate of an atom contributes an amount $Mg\delta z$ to the energy change $\Delta E$ of the system entering the MC acceptance probability $P_{acc}=min[1,exp(-\Delta 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 $\mu$m or larger, yielding $P_{c1}<0.003$ bar. The system size $L\simeq 315$ Å  used in our simulations, however, requires to include the linear regime to cover the pressure range smaller then $P_{c1}\simeq 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\eta)/(2-\eta)}+Bh^{3+2\mu}\simeq 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\rightarrow 0$ and $h\rightarrow\infty$:

[TABLE]

It has two fitting parameters $A$ and $B$ of which the latter is related to the elastic moduli by $B=(\kappa/k_{B}T)^{\mu}\tilde{g}Y/L^{4+2\mu}$.

The upper panel in Fig. 3 shows that for pressures up to $\sim$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}\simeq Bh^{3}$, which according to our analysis should apply for $P>P_{c3}$. This suggests that $P_{c3}\simeq$ 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\simeq 0.128$, implying a (bare) Young elastic modulus $Y=L^{4}B/g\simeq 299$ N/m. This is the value after the mentioned adjustment of $g_{\nu}$ 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 $\kappa=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 $\kappa$) can alternatively be used to determine $\tilde{g}$ 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\mu}=Y_{R}h^{3}/(gL^{4})$ to obtain the strain dependent $Y_{R}$:

[TABLE]

valid for $\epsilon<\epsilon_{c3}$, where $\epsilon_{c3}\simeq 0.005$ is the strain required to suppress anharmonicity, i.e. the strain beyond which the normal $P\sim h^{3}$ is applicable. Notice that this approach does not require knowledge of $\kappa$ nor $\tilde{g}$. For point load, the same approach can be used, but with an additional factor $\pi/(4L^{2})$ multiplying the right-hand side of eq. 13.

The right bottom panel shows a deviation from the $P\sim 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 Gao . The corresponding critical pressure, $P_{c4}$, should depend on the size as $P_{c4}=g_{4}/L$, with $g_{4}\simeq 25.2$ N/m derived from the value $P_{c4}\simeq 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, $\epsilon_{c4}=g_{\epsilon}h_{c4}^{2}/L^{2}$, the size dependence of $P_{c4}$ follows directly from eq. 4, neglecting the linear term, yielding $P_{c4}\simeq(gY/L)(\epsilon_{c4}/g_{\epsilon})^{3/2}$, with $\epsilon_{c4}\simeq 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 $\sim$0.2, a value indeed close to the strain where breaking was actually observed in our simulation, at a pressure $P_{br}\simeq 50$ kbar. This value of the breaking strain is similar to that found in a simulation study of graphene nanoribbons under uniaxial strain Lu . 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}\simeq$ 32 kbar for the simulated system size with a corresponding strain of $\sim$0.13, a graphene drum of 1 m in diameter gives a breaking force of $F_{c4}=\pi L^{2}P_{c4}/4\simeq 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 $\sim$ 0.328 due to the different values for $g$ and $g_{\epsilon}$, 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}/(\rho mg)\simeq 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 $\sigma_{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 $\pi\sigma_{0}h$. As this is an order $L^{2}$ larger then the linear contribution $~{}\kappa h/L^{2}$ arising from the FvK equations, its contribution increases with system size as $L^{2\eta-1}$, and can be significant as compared to the cubic term $~{}Yh^{3}/L^{2}$, even for $\mu$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 $\sigma_{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.

Acknowledgements.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 696656 – GrapheneCore1. We thank Cristina Gomez-Navarro, Julio Gomez-Herrero and Guillermo L/’opez-Pol/’in for interesting discussion.