Quantum elasticity of graphene: Thermal expansion coefficient and specific heat
I.S. Burmistrov, I.V. Gornyi, V.Yu. Kachorovskii, M.I. Katsnelson, A.D. Mirlin
aa r X i v : . [ c ond - m a t . m e s - h a ll ] N ov Quantum elasticity of graphene: Thermal expansion coefficient and specific heat
I. S. Burmistrov,
1, 2
I. V. Gornyi,
1, 3, 4, 5
V. Yu. Kachorovskii,
M. I. Katsnelson, and A. D. Mirlin
1, 3, 4, 7 L. D. Landau Institute for Theoretical Physics, Kosygina street 2, 119334 Moscow, Russia Laboratory for Condensed Matter Physics , National Research University Higher School of Economics, 101000 Moscow, Russia Institut f¨ur Nanotechnologie, Karlsruhe Institute of Technology, 76021 Karlsruhe, Germany Institut f¨ur Theorie der kondensierten Materie, Karlsruhe Institute of Technology, 76128 Karlsruhe, Germany A. F. Ioffe Physico-Technical Institute, 194021 St. Petersburg, Russia Radboud University, Institute for Molecules and Materials, NL-6525AJ Nijmegen, The Netherlands Petersburg Nuclear Physics Institute, 188300, St.Petersburg, Russia (Dated: October 2, 2018)We explore thermodynamics of a quantum membrane, with a particular application to suspendedgraphene membrane and with a particular focus on the thermal expansion coefficient. We showthat an interplay between quantum and classical anharmonicity-controlled fluctuations leads tounusual elastic properties of the membrane. The effect of quantum fluctuations is governed by thedimensionless coupling constant, g ≪
1, which vanishes in the classical limit ( ~ →
0) and is equalto ≃ .
05 for graphene. We demonstrate that the thermal expansion coefficient α T of the membraneis negative and remains nearly constant down to extremely low temperatures, T ∝ exp( − /g ).We also find that α T diverges in the classical limit: α T ∝ − ln(1 /g ) for g →
0. For grapheneparameters, we estimate the value of the thermal expansion coefficient as α T ≃ − .
23 eV − , whichapplies below the temperature T uv ∼ g κ ∼
500 K (where κ ∼ T ∼ − K. For
T < T , the thermal expansion coefficient slowly (logarithmically) approacheszero with decreasing temperature. This behavior is surprising since typically the thermal expansioncoefficient goes to zero as a power-law function. We discuss possible experimental consequencesof this anomaly. We also evaluate classical and quantum contributions to the specific heat of themembrane and investigate the behavior of the Gr¨uneisen parameter. PACS numbers: 72.80.Vp, 73.23.Ad, 73.63.Bd
I. INTRODUCTION
The thermal expansion coefficient α T is one of the mostimportant thermodynamic characteristics of any mate-rial. It is well known that α T is determined by an anhar-monicity of interatomic potentials binding atoms into acrystalline lattice [1–6]. Indeed, for a harmonic oscillator,the averaged displacement of the coordinate equals zeroindependently of the oscillation amplitude. By contrast,for an anharmonic oscillator, an averaged displacementdepends on the oscillation amplitude and, consequently,on temperature.For the most of materials, the thermal expansion co-efficients are positive, α T >
0. On the other hand, someexotic systems—including, in particular, complex metaloxides, polymers, and zeolites—are known to contractupon heating [7]. The interest to materials with α T < α T for all experi-mentally accessible temperatures: from very high (a fewhundreds of Kelvin) temperatures down to extremely lowtemperatures, T ∼ − K. Only at exponentially smalltemperatures,
T < T , the thermal expansion coefficient α T goes to zero. Since measurements of the elasticityof free-standing graphene have recently become accessi-ble to experimental techniques [22–24], the prediction ofan approximately temperature-independent negative α T can be verified experimentally. Our consideration is ap-plicable also to other two-dimensional materials.The fact that NTE is a natural property of layeredand 2D materials due to existence of bending (flexuralacoustic) mode which is anomalously sensitive to the de-formation has been recognized long time ago [25] anddicussed in a number of recent publications [19, 26–28]). In particular, graphite is known to have NTEat not too high temperatures [29], which may be ex-plained quantitatively within a quasiharmonic theory [2–6] via first-principles calculations of phonon spectra andGr¨uneisen parameters [26]. The same calculations pre-dicted that graphene should have NTE at any tempera-tures. More recently, atomistic simulations of graphene[30] confirmed the NTE at moderately high tempera-tures and also showed that the thermal expansion co-efficient changes sign with growing temperature due toanharmomic effects beyond the quasiharmonic approxi-mation. Qualitatively similar results were obtained inRefs. [31, 32]. Experiments confirm that the thermal ex-pansion coefficient of graphene at room temperature isnegative, with the absolute value as large as 0.1 eV − [33]. With the temperature increase up to 400 K , α T de-creases in absolute value [34], which may be consideredas a partial confirmation of the prediction [30]. At veryhigh temperatures, the thermal expansion coefficient ofgraphene is definitely positive but its precise measure-ment is difficult [35].The main focus of the present paper is on the range ofrelatively low temperatures. The behavior of the thermalexpansion coefficient of graphene (or, more generally, ofa 2D membrane) in this regime represents a challengingtheoretical problem. The microscopic Gr¨uneisen parame-ter for the bending mode is divergent at the phonon wavevector q → − /q (see Eq. 4 of [27]), which means that relevant phononsdetermining the thermal expansion are always classical(i.e, their energy is smaller than the temperature). Asa result, within the quasiharmonic approximation, thethermal expansion coefficient remains constant down toarbitrary low temperatures [27]. This is in stunning con-trast with the conventional behavior, α T → T → , which is usually associated with the third law of thermo-dynamics in view of the identity [36] (cid:18) ∂V∂T (cid:19) P = − (cid:18) ∂S∂P (cid:19) V (1)where V, P, T, S are volume, pressure, temperature, andentropy, respectively.Clearly, one may expect that quantum effects modifythe low- T behavior of α T . Quantum corrections to ther-modynamic properties have been calculated in Ref.37.Based on perturbative analysis, it was suggested thatclassical theory becomes inapplicable already at reason-ably high temperatures, about 70–90 K and that thethermal expansion coefficient goes to zero in a conven-tional way—i.e., as a power law—at T → . Here we willshow that the situation is, in fact, much more exotic andthe classical expression for the thermal expansion coef-ficient [27] remains valid, with some logarithmic correc-tions, till very low temperatures. The thermal expansioncoefficient does go to zero at zero temperature but slowerthan any power of the temperature.In this paper, we explore systematically thermody-namic properties of a graphene membrane, with a par-ticular focus on the thermal expansion coefficient α T , inthe whole range of temperatures. We show, in particular,that the behavior of α T in the limit T → g [definition of g is given by Eq. (12)below], which vanishes in the classical limit ~ → . As a one of the most important results of the paper, wedemonstrate that the thermal expansion coefficient of themembrane remains negative and nearly constant downto all realistic temperatures. For graphene parametersthis constant value is about − .
23 eV − . For a genericmembrane, we find that in the limit of weak coupling, thevalue of α T depends logarithmically on the bare value of the coupling g (the value of g at the atomic scales): α T ∝ − ln(1 /g ) , for g ≪ , and consequently diverges in the classical limit ~ → . We also demonstrate that, with decreasing temperature, α T starts to approach zero only when T drops below anexponentially small temperature scale T T ∼ g κ e − /g . (2)Here κ is the bare value of the bending rigidity of themembrane (equal to ≃ | α T | at such exponentially low temperaturesis logarithmically slow.We also calculate the specific heat. In particular, weshow that the leading contribution to both C V and C P isdetermined by classical effects, while the difference C P − C V is proportional to g and, thus, has a purely quantumnature.The outline of the paper is as follows. Section II con-tains a qualitative discussion of the physics of fluctuation-induced elasticity of a membrane. In Sec. III we presentthe classical and quantum renormalization-group (RG)formalism as well as basic equations for the observablesof interest. In Sec. IV we use the results of the precedingsection to calculate the temperature dependence of thethermal expansion coefficient and of the specific heat.Our results are summarized in Sec. V. Various technicalaspects of our analysis are presented in Appendices A–D. II. FLUCTUATION-INDUCED ELASTICITY OFA MEMBRANE
A free-standing graphene is a remarkable example ofa 2D crystalline membrane [19, 28]. Elastic propertiesof such a membrane are characterized by the bendingrigidity κ (which is quite high for graphene, κ ≃ Y and the bulk modulus B : Y = 4 µ ( µ + λ )2 µ + λ , B = µ + λ, (3)where µ and λ are the Lame coefficients (for graphene λ ≃ · ˚A − and µ ≃
10 eV · ˚A − , see Ref.30). A dis-tinct feature of a free-standing 2D membrane is the ex-istence of out-of-plane phonon modes—flexural phonons(FP) [19, 25, 28, 38]. In contrast to in-plane phononswith the conventional linear dispersion, the FP are verysoft, ω q ∝ q . As a consequence, the out-of-plane ther-mal fluctuations are unusually strong and tend to destroya membrane by driving it into the crumpled phase [38].The tendency of the membrane to crumpling is at theheart of the mechanism leading to NTE. Let us explainthis in more detail. The vector r describing point at themembrane surface depends on the 2 D coordinate x thatparametrizes the membrane and can be split into threeterms r = ξ x + u + h , (4) membrane with fluctuations membrane without fluctuations FIG. 1: Temperature-induced shrinking of a graphene mem-brane. Due to classical and quantum fluctuations, the mem-brane shrinks in the longitudinal direction (bottom), as com-pared to the membrane without fluctuations (top). where vectors u = u ( x , t ) , h = h ( x , t ) represent in-planeand out-of-plane phonon fields, respectively. The globalstretching factor ξ is equal to unity at zero temperaturebut gets reduced with increasing T due to out-of-planefluctuations (see Fig. 1), becoming zero at the crumplingtransition temperature T cr . Remarkably, the thermal ex-pansion of membrane is nonzero even if one neglects theanharmonicity of in-plane and out-of-plane modes as wellas anharmonic coupling between them [25]. Within suchan approximation the only relevant anharmonicity is dueto the coupling of global stretching and FP. Treating FPas classical random fields yields [39] ξ = 1 − h ∂ α h ∂ α h i / , (5)which results in a negative, logarithmically divergentvalue of α T = L − dA/dT : α T = ∂ξ ∂T = − π κ ln (cid:18) LL uv (cid:19) . (6)Here, A = ξ L is the sample area, L is the system size atzero temperature, and L uv is the ultraviolet cutoff of theorder of the lattice constant. Evidently, if this equationwere fully correct, the membrane would not exist in thethermodynamic limit, L → ∞ . Actually, the anharmoniccoupling between u and h fields leads to a renormaliza-tion of the bending rigidity, which becomes momentum-dependent, κ → κ q and flows to infinity, κ q → ∞ , for q → α T with the system size, thus yielding a finitevalue for the crumpling transition temperature, T cr < ∞ ,for L → ∞ . Below T cr one gets [39, 57] ξ = 1 − T /T cr , (7)and, consequently, a finite negative value of the thermalexpansion coefficient, α T = − /T cr . These arguments indicate that the role of anharmonic-ity in a free-standing 2D membrane is remarkably differ-ent as compared to the case of a 3D crystal (or to thecase of graphene on substrate). Indeed, in the 3D case, the anharmonic coupling between phonons determines anonzero value of α T (which is zero in the harmonic ap-proximation). In contrast, in a free-standing 2D mem-brane, such coupling leads to a suppression of an infinitevalue of α T predicted within the harmonic description ofin-plane and flexural modes down to a finite value.In order to find T cr and α T , one should investigate therenormalization from the ultraviolet energy scale downto the infrared scale. Such renormalization was inten-sively discussed more than two decades ago [38, 40–56]in the purely classical approximation (under assumptionthat the temperature is higher than the relevant frequen-cies of in-plane and flexural modes) in connection withbiological membranes, polymerized layers, and inorganicsurfaces. The interest to this topic has been renewedmore recently [37, 58–68] after discovery of graphene.It was found [40–47, 53] that the anharmonic couplingof in-plane and out-of-plane phonons leads to power-lawrenormalization of the bending rigidity, κ → κ q , where κ q ≃ κ (cid:18) q ∗ q (cid:19) η , for q ≪ q ∗ , (8)with a certain critical exponent η . Here the momentumscale (see e.g. Eq. (154) of Ref. [39]) q ∗ = r d c Y T π κ ∼ √ d c µT κ (9)is the inverse Ginzburg length which separates the re-gions of conventional ( q > q ∗ ) and fractal ( q < q ∗ ) scal-ing [19, 28, 38, 40]. In other words, the anharmonicityof flexural modes becomes important at q < q ∗ . Thelast expression in Eq. (9) is an order-of-magnitude es-timate where we have discarded numerical coefficientsas well as a difference in values of the elastic moduli, Y ∼ B ∼ µ ∼ λ . We will present such estimates also inseveral cases below in order to emphasize the scaling ofobservables with parameters of the problem.The critical exponent η was determined within severalapproximate analytical schemes [43, 45, 46, 53, 58]. Inparticular, for a 2D membrane embedded into a spaceof large dimensionality d ≫
1, one can find analytically η = 2 /d c ≪ , where d c = d − . Numerical simulations for physical 2D membrane embed-ded in 3D space ( d c = 1) yield η = 0 . ± .
10 [51] and η = 0 . ± .
04 [56]. Atomistic Monte Carlo simulationsfor graphene gives the value η ≈ .
85 [28, 65]; approx-imately the same value has been derived via functionalrenormalization group approach [58].In the limit d c ≫ , the critical temperature of thecrumpling transition can be also calculated analyticallywithin a classical approximation, yielding [39, 57] T cr = 4 πη κ d c = 8 π κ d . (10)As a consequence, the thermal expansion coefficient isnegative and independent of temperature α d c →∞ T = − d c πη κ . (11)Evidently, one expects that this result should fail at lowenough temperatures due to the quantum effects.Some aspects of this problem were discussed recently[27, 37, 63, 67]. One scenario of the suppression of α T is the emergence of a “mass” (term quadratic in momen-tum q ) in the propagator of flexural phonons or, equiv-alently, 1 /q divergence of the effective bending rigidity[37]. This mass arises naturally within perturbative cal-culations [37]. However, as it will be shown below by anexplicit calculation of the free energy, it is nothing else asthe full tension (as was pointed out in Ref. 63) and thusis zero for a free membrane (without external stress).As was shown in Ref. 62, quantum fluctuations mayalso lead to renozmalization of elastic coefficients. Thisrenormalization can be described [62, 64] in terms of aflow of the dimensionless quantum coupling constant g = 3(6 + d c )128 π ~ Yρ / κ / ∼ ~ µρ / κ / , (12)where ρ is the mass density of membrane. Since g ∝ ~ , it vanishes in the classical limit ~ → . For graphene,the bare value g of this constant at the ultraviolet scale, q ∼ q uv , on the order of the lattice constant is quite small, g ≃ /
20. Physically, this happens because g contains ρ and, consequently, the atomic mass in the denominator.In the sequel, we develop a theory of thermal expan-sion of an elastic membrane that includes both classicaland quantum effects. We show that there is an additionalcontribution to α T , which originates from the region ofmomenta q ∗ < q < q ∗ / √ g and is not taken into accountin Eq. (11). This contribution is logarithmically large[ ∝ − ln(1 /g )] for small coupling g. Quantum fluctuationsoriginate from q > q ∗ / √ g. Evidently, such interval of q exists only when q ∗ / √ g < q uv . The temperature foundfrom the condition q ∗ / √ g ∼ q uv is given by T uv ∼ g κ ( ∼
500 K for graphene). This temperature is determinedby the bare value g of the quantum coupling constantand plays the role of the Debye temperature. For T < T uv the problem becomes quantum in terms of statistics, i.e.,some phonons are frozen out. On the other hand, theeffect of quantum fluctuations (i.e., those with momentain the range q ∗ / √ g < q < q uv ) remains negligibly smallin a wide temperature interval T < T < T uv , where T ∼ T uv exp( − /g ). Therefore, from the point of viewof fluctuation-induced renormalization, the problem re-mains classical down to the temperature T uv , which isexponentially small for g ≪
1. Only at very low tem-peratures,
T < T , quantum fluctuations come into playand, as a result, the thermal expansion coefficient getslogarithmically suppressed.Our goal in this paper is to develop a theory of thermalexpansion valid for all temperatures in the range T < T uv where the main contribution to α T originates from theflexural phonons and therefore is negative. This requiresa development of formalism incorporating both classicaland quantum renormalization effects, which is the sub-ject of the next Section. III. FORMALISMA. Thermodynamics of an elastic membrane
We consider a generic 2 D membrane embedded in the d -dimensional space ( d > L ( { r } ) = ρ ˙ r + κ r ) (13)+ µ ∂ α r ∂ β r − δ αβ ) + λ ∂ γ r ∂ γ r − , which can be obtained from the general gradient expan-sion of elastic energy [42] by using a certain rescaling ofcoordinates (see discussion in Ref. 39). The subscript0 in notations for elastic coefficients in Eq. (13) meansthat these are bare values at q ≃ q uv , where q uv is theultraviolet cut off. The d − dimensional vector r = r ( x , τ )is given by Eq. (4) with u = ( u , u ) , h = ( h , ..., h d c ),while ˙ r = d r /dτ, where τ is imaginary time.The strategy of calculations is as follows. We assumethat the ahrahmonic phonon interaction leads to a renor-malization of elastic coefficients κ → κ q , µ → µ q , λ → λ q . Hence, in order to calculate the free energy,we replace Eq. (13) with a harmonic Lagrangian den-sity containing renormalized elastic moduli. The detailsof calculation are relegated to Appendice A and B. Theobtained free energy has the form [see Eq. (B14)] FL = − σ B + σ ( ξ −
1) (14)+ d c X q ω ln (cid:0) κ q q + σq + ρω (cid:1) , where σ = 1 L ∂F∂ξ , (15)is the external stress, P q ω stands for T P ω R d q / (2 π ) , andthe summation P ω runs over bosonic Matsubara frequen-cies. This approximation is in spirit of self-consistentphonon theory [69] where, in analogy to Landau Fermiliquid theory for fermions, it is supposed that the en-tropy is renormalized by phonon-phonon interaction onlyvia the change of their dispersion relation and phonondamping is neglected. One can assume that, within somenumerical factors, this gives correct temperature depen-dences of all thermodynamic quantities.Since FP are much softer than in-plane modes, we ne-glected contribution of in-plane modes in Eq. (14). Al-though this approximation looks quite natural, it needssome justification. Indeed, as was found in Ref. 37, theanharmonicity induces a small (in adiabatic parameter)ultraviolet-divergent contribution σ (“built-in tension”)to the coefficient in front of the q term in the propaga-tor of FP. Such term would suppress α T (in a power-lawway) at low temperatures. In fact, this term is exactlycancelled by another contribution, which was overlookedin Ref. 37. Technically, this additional contribution arisesdue to the coupling between in-plane phonons and globalstretching. Here, we discuss the problem on the qualita-tive level, relegating the details of calculations to Appen-dices A and B. Substituting Eq. (4) into strain tensor ofthe membrane u αβ = 12 ( ∂ α r ∂ β r − δ αβ ) , and leaving only linear with respect to fluctuation terms,we find that, in the harmonic approximation, the straintensor is proportional to the global stretching and to thegradient of in-plane deformations: u harmonic αβ = ξ ( ∂ α u β + ∂ β u α ) / . (16)The energy of the in-plane fluctuations E in − plane isquadratic with respect to u αβ [38] and, as follows fromEq. (16), is proportional to ξ . Hence, there exists ancontribution to the stress δσ = h E in − plane i ξ L , (17)where averaging is taken with the action correspondingto the Lagrangian (13). As shown in Appendixes A andB δσ = σ . (18)This implies that the quadratic-in- q part of the self-energy of FP is given by ( σ + σ − δσ ) q = σq , where σ is the external tension. In other words, coupling of in-plane modes to the global stretching ξ leads to additionalcontribution to the FP’s self-energy which exactly can-cels the quadratic correction arising due to anharmoniccoupling of in-plane modes with FP.This cancelation has a deep physical meaning. Thereare two different definitions of the tension σ . First, it canbe obtained from the standard thermodynamic relation,as a derivative of the free energy with respect to systemvolume, see Eq. (15). Second, the tension can be foundas a coefficient in front of the quadratic-in- q term in theFP propagator, σ = ∂G − q ∂ ( q ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) q =0 . (19)The equivalence of two definitions, while quite transpar-ent physically, represents a very non-trivial Ward iden-tity [44, 47]. We explicitly verify this identity within theone-loop RG analysis in Appendixes A and D. FIG. 2: Characteristic momentum scales in the problem. Re-gions of quantum and classical RG are shown, along with thecorresponding RG equations for bending rigidity. In theseequations, η and θ are the classical and quantum critical in-dices, respectively, while g is the quantum coupling constant. In-plane modes lead also to a small ultraviolet renor-malization of ξ which we neglect here (see Appendix Bfor more detail).The relation between ξ and σ is found from the con-dition ∂F/∂σ = 0 , which yields, after summation overMatsubara frequencies, ξ − ξ = σB + s ( σ ) , (20)where ξ = 1 − d c πρ q uv Z dq q ω q coth ω q T ! , (21)is the longitudinal deformation of membrane in the ab-sence of stress, while s = d c πρ q uv Z dqq " coth (cid:0) ω q / T (cid:1) ω q − coth ( ω q / T ) ω q , (22)is the stress-induced correction ( s σ → = 0), which leads,in particular, to anomalous Hooke’s law [57]. Here ω q and ω q are FP frequencies for stressed and unstressedmembrane, respectively, ω q = s κ q q + σq ρ , ω q = r κ q ρ q . (23)In order to complete the calculation of α T , we shouldfind renormalization of κ in the whole interval of mo-menta, 0 < q < q uv , and substitute the renormalizedfunction κ q into Eqs. (20)-(23). This will allow us to de-termine the global stretching factor as a function of theapplied stress and temperature, ξ = ξ ( σ, T ) , and thus toevaluate the thermal expansion coefficient, α T = (cid:20) ∂ξ ( σ, T ) ∂T (cid:21) σ . (24)Hence, we focus in the next subsection on renormaliza-tion of elastic constants. B. Renormalization group
In this section, we find the flow of the elastic moduliwith decreasing q down from q uv by using the pertur-bative RG approach (for a recent discussion of a quan-tum non-perturbative RG see Ref. 68). The classical andquantum RG are separated by q = q T found from thecondition ~ ω q ≃ T : q T ∼ q ∗ √ g . (25)For g ≪ , the characteristic scales of the problem areshown in Fig. 2. At large spatial scales (for q ≪ q ∗ ),a classical RG apply, so that the bending rigidity scalesaccording to Eq. (8). In the interval q ∗ < q < q T , the FPfrequency is still small compared to T , so that classicalapproach remains applicable. However, the renormaliza-tion of κ in this region is small [38], ( κ q − κ ) / κ ∼ q ∗ /q .The quantum RG operates in the interval between q T and q uv , with the former scale serving as an infrared and thelatter as an ultraviolet cutoff. In order to derive quan-tum RG equations, we notice that for q < p µ / κ , theFP are much softer than the in-plane modes: ω q ≪ ω ⊥ , k q , (26)where ω ⊥ q = q p µ/ρ, ω k q = q p (2 µ + λ ) /ρ are frequencies of the transverse and longitudinal in-plane phonos respectively. We further notice, that forgraphene, the value of p µ / κ is on the order of the in-verse lattice constant. This implies that one can use thisvalue as the ultraviolet atomic momentum scale: q uv ∼ p µ / κ . (27)In view of Eq. (27), the interval q T < q < q uv for quan-tum RG exists provided that the temperature is not toohigh, T < T uv , where T uv ∼ g κ . (28)For higher temperatures one can fully neglect quantumeffects. Equation (28) implies that for graphene ( g ≃ / κ ≃ T uv is of the orderof 500 K. This temperature plays the role of the Debyetemperature in our model.As follows from Eq. (26), the retardation effects canbe neglected and the quantum RG flow can be found ina full analogy with classical RG, [39] where FP field isconsidered to be static. Technical details are describedin Appendices C and D where two alternative derivationsof RG equations are presented. For in-plane moduli, onegets the following RG equations: dd Λ 1 Y = 38 dd Λ 1 B = d c g (6 + d c ) Y , (29) where Λ = ln ( q uv /q )is the logarithm of the running RG scale q , and the cou-pling g is given by Eq. (12). For d c = 1 , these equationsare equivalent to Eqs. (9) and (10) of Ref. 62. As seen,there is invariant subspace of elastic moduli, Y = 8 B/ , which is conserved by the RG flow.The RG flow of the bending rigidity κ is given by equa-tion d κ d Λ = 4 g κ d c + 6 , (30)which agrees up to a sign with Eq. (11) of Ref. 62. [Thesign was recently corrected in Ref. 64.] From Eqs. (12),(29), and (30), we find dgd Λ = − g , (31)which again agrees with the result of Ref. 62 once thesign error is corrected [64]. (We notice that the numer-ical coefficient in definition of g in Ref. 62 is different).The negative sign in Eq. (31) is of key importance: it im-plies that the quantum anharmonicity effects stabilize themembrane increasing the effective bending rigidity. Inother words, the flat phase of the membrane is perfectlydefined in the limit of zero temperature and infinite sys-tem size. This is crucially important for low-temperaturebehavior of the α T It is the growing effective bendingrigidity at low temperatures which suppresses the ther-mal expansion coefficient. However, contrary to all pre-viously known situations this provided only logarithm-in-power vanishing of the thermal expansion coefficientat T → g , the bending rigidity κ ,and the in-plane moduli Y and B : g = g g Λ , κ = κ (1 + g Λ) θ , (32) Y = Y (1 + g Λ) − θ/ , B = 1 B + 83 (cid:18) Y − Y (cid:19) . (33)Here θ = 4 d c + 6 (34)is the quantum anomalous exponent. For graphene (or,more generally, for a 2D membrane in a 3D space), d c = 1and θ = 4 / . The bare coupling constant g is quite smallfor graphene (about 1/20). With increasing spatial scale,the running coupling constant g decreases according toRG equation (32) from g (equal to ≃ g ≪ . This justifies[62] the applicability of one-loop RG approach.The quantum RG stops at q = q T . The overall pictureof the renormalization of κ is as follows. The RG flowstarts at q ∼ q uv , where κ has a bare value κ . In the in-terval q T < q < q uv , the bending rigidity grows accordingto Eq. (32). The value of κ q at the edge of the quantuminterval (at q ≃ q T ) and the corresponding value of g aregiven by κ = κ (Λ T ) , g = g (Λ T ) , where κ (Λ) and g = g (Λ) are given by Eq. (32) andΛ T = ln ( q uv /q T ) ≈ ln p T uv /T. Below the couplings without indication of momentumscale (such as g , κ , B , Y ) will be understood as de-fined on the scale Λ T governed by the temperature. It isworth noting that values of q ∗ and q T are determined bythe renormalized elastic moduli: q T ∼ √ Y T κ √ g ∼ √ Y T κ √ g (1 + g Λ T ) θ/ ; (35) q ∗ ∼ √ Y T κ ∼ √ Y T κ g Λ T ) (2+ θ ) / . (36)In the interval q ∗ < q < q T , the bending rigidity does notchange essentially (i.e., it changes by a factor of orderunity). Finally, at lowest momenta q < q ∗ , the bendingrigidity scales according to Eq. (8) with κ = κ (Λ T ) . Now we are in a position to calculate the integrals en-tering Eqs. (21) and (22) and to get final formulas gov-erning the thermodynamics of graphene. Results of thisanalysis are presented in the next section.
IV. RESULTSA. Thermal expansion coefficient: Zero tension
The main contribution to ξ as given by Eq.(21) comesfrom the region q < q T . For η ≪ , a simple calculationyields ξ ≈ − d c T π κ (cid:20) η + ln (cid:18) g (cid:19)(cid:21) . (37)Here we neglect terms of the order of g as well as onesof the order of T / κ coming from q > q T . We, thus, ob-tain the thermal expansion coefficient of an unstressedmembrane: α T ≈ − d c π κ (cid:20) η + ln (cid:18) g (cid:19)(cid:21) . (38)Two terms in the square brackets represent contributionsof momentum intervals q < q ∗ and q ∗ < q < q T , respec-tively. Comparing Eq. (38) with Eq. (11), we observe twodifferences.Firstly, a logarithmic-in- g term (reflecting the contri-bution of the momenta q ∗ < q < q T ) appeared in thesquare brackets. This term can be neglected for a generic membrane embedded in the space of high dimensionality( d c → ∞ , η →
0) [57]. However, for graphene, where η ≃ . /g ≃
20, the two terms give comparablecontributions. On the other hand, for a “nearly classi-cal” 2D membrane in 3D space that has the same η andmuch larger 1 /g , the logarithmic contribution would bedominant.Secondly, α T becomes now a slow function of temper-ature. This dependence deserves a special attention. Asfollows from Eqs. (32) and (38), α T remains negative andnearly constant in an extremely wide temperature range: α T ≈ α max = − d c π κ (cid:20) η + ln (cid:18) g (cid:19)(cid:21) , T ≪ T ≪ T uv . (39)Here T is given by Eq. (2), yielding for graphene T ≃ − K. Thus the thermal expansion coefficient ofgraphene remains nearly constant within almost twentydecades of temperature! Using graphene parameters( d c = 1 , κ ≃ α T ≃ − .
23 eV − .Only at exponentially low temperatures, T ≪ T , thethermal expansion coefficient starts to decay logarithmi-cally with decreasing temperature: α T ≃ − d c π κ ln ln( T uv /T )[( g /
2) ln( T uv /T )] θ , T ≪ T . (40)The temperature dependence of α T is shown in Fig. 3 for T ≪ T uv .It is instructive to analyze how the classical limit( ~ →
0, implying g → α T depends logarithmically on g . Conse-quently, the thermal expansion coefficient (39) divergesin the classical limit. It should be emphasized, however,that the range of validity of Eq. (39) shrinks in this limitin view of Eq. (28). For temperatures above T uv , thelogarithmic term in Eq. (39) gets modified, becomingtemperature-dependent: α T ≈ − d c π κ (cid:20) η + ln (cid:18) q uv q ∗ (cid:19)(cid:21) , T ≫ T uv . (41)In other words, the function | α T ( T ) | has a maximumat T ≃ T uv (this maximum is not shown in Fig. 3,which is plotted for T ≪ T uv ). That is, it is themaximal value α max that diverges in the classical limit: α max ∝ − ln(1 /g ) . B. Thermal expansion coefficient: Finite tension
We turn now to a generalization of the above resultsto the case of non-zero tension, σ = 0 . In this case, thereappears a new characteristic momentum q σ determinedby the condition κ q q σ ∼ σ. FIG. 3: Temperature dependence of the thermal expansioncoefficient for unstressed ( σ = 0) and stressed ( σ = 0) mem-brane. The tension σ suppresses the thermal expansion at T < T σ . This momentum increases with σ , reaching the value of q ∗ for σ = σ ∗ , where σ ∗ ∼ κ q ∗ q ∗ ∼ d c Y T κ . (42)We remind the reader that, for temperatures below T uv ,couplings without the momentum indicated [such as Y and κ in Eq. (42)] are understood as those includingthe quantum renormalization, i.e., defined on the tem-perature scale q T . We have also taken into account inEq. (42) that there is no essential renormalization of κ between the scales q T and q ∗ .The physical meaning of σ ∗ was discussed in Ref. 57.For σ > σ ∗ the membrane shows linear Hooke’s law, [57],while for σ < σ ∗ the stress-strain relation is of a power-law form, ξ − ξ ∝ σ α , with an anomalous exponent α given by α = η − η . One can express σ ∗ in terms of the bulk modulus σ ∗ = d c CBT π κ , (43)where C is a numerical coefficient of order unity, which ischosen from the requirement that the low-stress deforma-tion takes the form given by upper line of Eq. (48) below.This coefficient is determined by relations between val-ues of elastic constants. (In principle, C depends on theratio of elastic module and, therefore, is a slow functionof temperature [70].) The room-temperature value of C for graphene, C ≃ .
25, was found in Ref. 57 from a com-parison of the analytic strain-stress relation with results of atomistic simulations of Ref. 71. [The definition ofthe coefficient C in Eq. (43) differs from that in Ref. 57by an additional factor 2 /η, which is ≃ σ up to the value σ ∗ /g, themomentum q σ reaches the boundary of the quantum re-gion: q σ ≃ q T for σ ≃ σ ∗ /g. (44)Since we want to analyze a temperature dependence ofmembrane properties at non-zero tension σ , it is usefulto introduce a characteristic temperature T σ determinedby the condition σ = σ ∗ ( T ), T σ = 4 π κ σd c CB ∼ κ σB . (45)The stress-induced deformation s ( σ ) can be separatedinto classical and quantum parts: s ( σ ) = s cl ( σ ) + s quant ( σ ) , where s cl ≈ d c T πρ q T Z dqq (cid:18) ω q ) − ω q (cid:19) , (46) s quant ≈ d c πρ q uv Z q T dqq (cid:18) ω q − ω q (cid:19) . (47)For σ < σ ∗ /g , the classical contribution is essentiallynon-perturbative with respect to σ (see Ref. 57): s cl ≈ σ ∗ B ( (1 /α ) ( σ/σ ∗ ) α , for σ < σ ∗ (1 / C )[2 /η + ln( σ/σ ∗ )] , for σ ∗ < σ < σ ∗ /g. (48)In contrast, the quantum contribution can be calculatedperturbatively by expansion to the leading order in σ : s quant = − d c σ πρ q uv Z q T dqq ∂ω − q ∂σ ! σ =0 = σ Λ T Z d Λ ∂B − ∂ Λ = σB − σB . (49)Here we have used Eq. (29). Substituting s quant intoEq. (20), we see that quantum effects lead to a simplerenormalization of the first term in the r.h.s. of Eq. (20): σ/B → σ/B. For a sufficiently large stress (or for a sufficiently lowtemperature), σ > σ ∗ /g , the momentum q σ becomeslarger than q T , so that one can neglect the term κ q q in Eq. (23) in comparison with σq in the whole classicalregion. In this case, Eq. (20) takes the form ξ = ξ σ + σB σ , where B σ is the renormalized value of the bulk modulus[see Eq. (33)] at the scale q σ ≃ p σ/ κ , and ξ σ ≈ − d c π √ ρσ Z dq q (cid:20) coth (cid:18) √ σq √ ρT (cid:19) − (cid:21) . Here we neglected small ( ∼ g ) terms. The integral in thisformula scales with decreasing temperature as T . Summarizing the obtained results, we find ξ − ≈ − d c T π κ (cid:20) η + ln (cid:18) g (cid:19)(cid:21) + σ ∗ αB (cid:18) σσ ∗ (cid:19) α , for σ < σ ∗ ,σB − d c T π κ ln (cid:18) σ ∗ σg (cid:19) , for σ ∗ < σ < σ ∗ /g,σB σ − d c ρT ζ (3)2 πσ , for σ ∗ /g < σ. (50)From Eqs. (24) and (50), we obtain the thermal expan-sion coefficient α T as a function of applied stress: α T ≈ − d c π κ /η + ln(1 /g ) − C ( σ/σ ∗ ) α , for σ < σ ∗ , ln( σ ∗ /σg ) , for σ ∗ < σ < σ ∗ /g, ζ (3) ρT κ /σ , for σ ∗ /g < σ, (51)where C = 4 C (1 − η ) /η and ζ (3) ≃ .
202 is the Reimannzeta function.The dependence of α T on T for stressed and unstressedmembrane following from Eq. (51) is shown in Fig. 3. (Inthis figure, we assume for the stressed case that g T σ ≫ T . ) At high temperatures, T ≫ T σ , the external tensionresults in a power-law dependence of α T on temperature: α T ≈ − d c π κ (cid:20) η + ln 1 g − C (cid:18) T σ T (cid:19) α (cid:21) . (52)Below T σ , the absolute value of the thermal expansioncoefficient decreases logarithmically: α T ≈ − d c π κ ln TgT σ , gT σ ≪ T ≪ T σ . (53)Finally, at still lower temperatures, T ≪ gT σ , we find α T ∝ T , see the third line in Eq. (51). C. Effective bulk modulus
Differentiating Eq. (50) over σ, we find expression foreffective bulk modulus1 B eff = (cid:18) ∂ξ ∂σ (cid:19) T (54)= (1 /B )( σ ∗ /σ ) − α , for σ < σ ∗ , /B + d c T / π κ σ, for σ ∗ < σ < σ ∗ /g, /B σ + d c ζ (3) ρT /πσ , for σ ∗ /g < σ, FIG. 4: Temperature dependences of (a) strain ξ for fixedstress σ , and (b) stress for fixed strain. Positions of crumpling(in the left panel) and buckling (in the right panel, for ξ < Two upper lines of Eq. (54) were obtained previously inRef. 57. In the second and third line we took into accountsmall temperature-dependent corrections. The third lineshows that B eff slowly increases with σ at σ > σ ∗ /g dueto suppression of quantum RG by external stress: B eff ≃ B σ = B [1 + ( g /
2) ln( µ /σ )] − θ/ . (55)[Here we assume for simplicity that Y = 8 B/ , useEq. (27) and write ln( q uv /q σ ) with the logarithmic pre-cision]. D. General phase diagram of membrane
In Fig. 4, we plot the temperature dependence of thestrain ξ for fixed stress σ (Fig. 4a) as well as the temper-ature dependence of the stress for fixed strain (Fig. 4b).0As seen from Fig. 4a, with increasing temperature for thefixed σ, the membrane undergoes crumpling transition( ξ → T cr ( σ ) canbe found from the upper line of Eq. (50) by requirement ξ = 0 and by taking into account Eqs. (43) and (32).The plot shown in Fig. 4b corresponds to an experimen-tal setup in which the membrane in-plane area is keptfixed, while the temperature is varied. If this area issmaller than the intrinsic zero-temperature area of themembrane (i.e., ξ < E. Specific heat
We evaluate now the temperature dependence of thespecific heat. Similar to the thermal expansion coefficient α T , both constant-volume ( C V ) and constant-pressure( C P ) specific heat capacities are determined by FP. Toevaluate them, we first determine the entropy of themembrane, S = − L − ( ∂F/∂T ) ξ , where F is the free energy given by Eq. (14). We find S = d c Z d q (2 π ) [( N q + 1) ln( N q + 1) − N q ln N q − ( N q + 1 / ∂ω q /∂T ) σ ] , (56)where N q = [exp( ω q /T ) − − is the Bose function. Itis worth noting that the last term in the square brack-ets in Eq. (56) is nonzero because of the temperature-dependent renormalization of the bending rigidity κ [seeEq. (8)]. By using Eq. (56), one can calculate the specificheat capacities, C P,V = − T ( ∂S/∂T ) P,V . A straightforward calculation shows that, in the leadingorder, both C P and C V are given by the classical formulafor phonons with a parabolic spectrum, C P ≈ C V ≈ d c Z d q (2 π ) ω q ∂N q ∂T ≈ πd c √ ρT √ κ . (57)In the presence of a finite tension, this result is validfor temperatures T ≫ gT σ (or, equivalently, the tension σ ≪ σ ∗ /g ). At low temperatures, T ≪ gT σ , the specificheat is proportional to T : C P ≈ C V ≈ ζ (3) d c ρπσ T . (58)While the leading contribution to both C P and C V is dueto classical fluctuations, the difference C P − C V is small,proportional to g, and, therefore, is due to the quantum effects: C P − C V C P = T ( ∂ξ /∂T ) σ C P ( ∂ξ /∂σ ) T (59) ∼ g ln (1 /g ) ( σ/σ ∗ ) − α , for σ ≪ σ ∗ , [ln( σ ∗ /gσ )] , for σ ∗ ≪ σ ≪ σ ∗ /g, ( B σ /B )( σ ∗ /gσ ) , for σ ∗ /g ≪ σ Hence, at high temperatures, T ≫ T σ , we find C P − C V C P ∼ g ln (cid:18) g (cid:19) (cid:18) T σ T (cid:19) − α . (60)At intermediate temperatures, gT σ ≪ T ≪ T σ , the resultreads C P − C V C P ∼ g ln TgT σ , (61)while at very law temperature T ≪ gT σ , we obtain C P − C V C P ∼ g (cid:18) TgT σ (cid:19) B σ B ∼ ρT B σ σ . (62)We emphasize that C P − C V vanishes in the absence ofthe external tension ( σ = 0), as follows from Eq. (60)with T σ = 0. F. Gr¨uneisen parameter
Finally, we consider the macroscopic Gr¨uneisen param-eter γ = α T ( ∂σ/∂ξ ) T C V , (63)which is an important characteristics of thermomechan-ical properties of a system. We find that the Gr¨uneisenparameter is negative for all values of stress: γ ∼ − g κ T ln(1 /g ) ( σ/σ ∗ ) − α , for σ ≪ σ ∗ , ln( σ ∗ /gσ ) , for σ ∗ ≪ σ ≪ σ ∗ /g, ( B σ /B )( σ ∗ /gσ ) , for σ ∗ /g ≪ σ. (64)Several points deserves special attention. First of all, wesee that the absolute value of γ has a maximum as afunction of σ at σ ≃ σ ∗ . On the other hand, for fixed σ, | γ | is a monotonously decreasing function of T. Mostimportantly, in the limit of low temperature, γ turns outto be non-zero − γ ∼ B σ σ , for T → . (65)Schematic dependence of the Gr¨uneisen parameter ontension and temperature is illustrated in Fig. 5.1 FIG. 5: Schematic dependence of the Gr¨uneisen parameteron tension at fixed temperature (a) and on temperature forfixed tension (b).
G. Third law of thermodynamics
Let us point out that the third law of thermodynam-ics manifests itself in this problem in a somewhat curiousway. Indeed, the entropy (56) vanishes in the limit T → σ even if the quantum renormalization effects areneglected. On the other hand, it is easy to see from thefirst line of Eq. (51) that α T | σ =0 remains finite in thelimit T → . At first glance, this may seem to contradictto Eq. (1) that yields α T = ∂S/∂σ. This apparent con-tradiction is resolved by noticing that non-analyticity of S at T = 0 , σ = 0 leads to non-commutativity of limits T → σ → α T ( T, σ ) . Quantumrenormalization effects restore the vanishing of the ther-mal expansion coefficient at T = 0 independently of theorder of the limits, which is conventionally considered asa manifestation of the third law. V. CONCLUSION
To summarize, we have developed a theory of ther-momechanical properties of a suspended graphene mem-brane. We have shown that at zero tension the ther- mal expansion coefficient α T of free-standing graphene isnegative and temperature-independent in a very broadtemperature range, see Eq. (39). The underlying physicsof the negative expansion is the global shrinking of thegraphene membrane in the longitudinal direction due toclassical transverse fluctuations.The second term in Eq. (39) for α T is governed by thedimensionless quantum coupling constant, g ≪
1. Thiscoupling constant vanishes in the classical limit ~ → α T ) and is equal to ≃ . g ensures that α T remains T -independent down to extremely low temper-ature T , Eq. (2). For graphene parameters, we esti-mate the value of the thermal expansion coefficient as α T ≃ − .
23 eV − , which applies below the temperature T uv ∼ g κ ∼
500 K (where κ ∼ T ∼ − K. For
T < T , the abso-lute value of the thermal expansion coefficient starts todecrease logarithmically slowly with decreasing temper-ature, Eq. (40), since quantum effects lead to increaseof the bending rigidity. Our results imply that, contraryto naive expectations, quantum fluctuations do not leadto to the melting or crumpling of a 2D crystal but in-stead stabilize the membrane due to enhanced role of theanharmonicity.A finite tension σ suppresses the thermal expansion at T < T σ , where T σ ∝ σ is the characteristic temperaturewhich separates regimes of conventional ( T < T σ ) andanomalous ( T > T σ ) elasticity, see Fig. 3.We have also evaluated the temperature dependence oftension in a graphene membrane placed into a frame ofa fixed size ξL . With lowering temperature, a membranewith ξ < C P and C V are dominated by classical effects andare given by a standard expression for phonons withparabolic spectrum. On the other hand, a small differ-ence C P − C V is due to quantum fluctuations and shows avery non-trivial behavior as a function of the ratio T /T σ ,see Eqs. (60), (61), and (62). The same ratio T /T σ de-termines the temperature and stress dependence of theGr¨uneisen parameter, which turns out to be negative forall temperatures and tensions, being monotonous func-tion of T (for fixed σ ) and showing a minimum as a func-tion of σ (for fixed T ) for σ ≃ σ ∗ . Our results demonstrate that 2D materials are dra-matically different from 3D ones where, according toGr¨uneisen law, the thermal expansion coefficient is pro-portoinal to the heat capacity and goes to zero as T → α T experimentally. In particular, the thermalexpansion coefficient can be measured from the temper-ature shift of Raman spectra [34]. An alternative (andpossibly an easier) way of determining α T is providedby studies of van der Waals heterostructures such as2graphene/hBN, graphene/MoS , etc. [72, 73]. In thissetting, large thermal expansions of 2D materials at lowtemperatures would result in a strong temperature de-pendence of lattice mismatch which can be seen via re-construction of Moire patterns [74–76] or via character-istics of graphene bubbles on a substrate [77]. VI. ACKNOWLEDGEMENT
We thank E. I. Kats, V. V. Lebedev, and K. S.Novoselov for useful discussions. The work was sup-ported by the Russian Science Foundation (grant No.14-42-00044). MIK acknowledges a support by NWO viaSpinoza Prize.
Appendix A: Calculation of free energy
In this Appendix we calculate the free energy corre-sponding to the Lagrangian (13). The partition functionreads Z = Z { D r } e − S , S = Z β dτ Z d x L , (A1)where β = 1 /T. In terms of u , h , and ξ introduced inEq. (4), the Lagrangian density becomes L = ρ u + ˙ h )+ µ + λ "(cid:18) ξ − K (cid:19) − K + L , (A2)where ˙ u = ∂ u /∂τ , ˙ h = ∂ h /∂τ,K = 2 h u αα i x ,τ = h ∂ α h ∂ α h + ∂ α u ∂ α u i x ,τ , (A3) u αβ = (cid:2) ξ ( ∂ α u β + ∂ β u α ) + ∂ α h ∂ β h + ∂ α u ∂ β u (cid:3) / , (A4) h· · · i x ,τ denotes the spatial and time averaging, h· · · i x ,τ = Z d x L Z β dτβ · · · , and L = κ h ) + (∆ u ) ] + µu ij + λ u ii . (A5)Equation (A4) coincides with the conventional expressionfor the strain tensor of the membrane provided that ξ = 1and the term ∂ α u ∂ β u in the square brackets (which is ofthe second order in ∂ α u and thus much smaller than thefirst term) is neglected. Within such an approximation,and neglecting also a small term κ (∆ u ) / L ( u , h ) coincides with the textbook ex-pression for elastic energy of a nearly flat membrane. [38]The term ( µ + λ ) K / h , u , and h u couplings) with zero transferred momentum and energy, q = 0 and Ω = 0 (zeromode). After combining this term with the analogousquartic zero-mode term coming from L , the zero-modecontribution can be safely neglected because it gives anegligibly small correction to the self-energies of bothflexural and in-plain phonons (see detailed discussion inthe Supplementary Material of Ref. 57)In the ( q , ω ) representation, the action we are left withreads S = Z q ω (cid:20) ρω + κ q | h q ,ω | + | u q ,ω | ) (A6)+ µ | u q ,ωαβ | + λ | u q ,ωαα | (cid:21) + µ + λ L T (cid:18) ξ − K (cid:19) . Here R q ω stands for T P ω R d q / (2 π ) and summation goesover Matsubara frequencies ω = 2 πnT. The next step isto integrate exp( − S ) over u and h . To carry out thisintegration, we first perform a Hubbard-Stratonovich de-coupling of the last term in Eq. (A7) by an integral overan auxiliary field χ ,exp " − L ( µ + λ )2 T (cid:18) ξ − K (cid:19) = L p π ( µ + λ ) T × Z dχ exp (cid:26) − (cid:20) χ µ + λ ) − iχ (cid:18) ξ − K (cid:19)(cid:21) L T (cid:27) . This yields e − S = L p π ( µ + λ ) T Z dχe − S χ (A7) × exp (cid:26) L T (cid:20) − χ µ + λ ) + iχ (cid:0) ξ − (cid:1)(cid:21)(cid:27) , where S χ = Z q ω (cid:20) ρω + κ q − iχ | h q ,ω | + | u q ,ω | )+ µ | u q ,ωαβ | + λ | u q ,ωαα | (cid:21) . (A8)We will first discuss what happens in the harmonic ap-proximation and later include the anharmonic couplingbetween the in-plane and out-of-plane modes. In the har-monic approximation, S χ simplifies, S χ = Z q ω (cid:20) ρω + κ q − iχ | h q ,ω | + | u q ,ω | )+ µξ q | u q ,ω | + ( µ + λ ) ξ | qu q ,ω | (cid:21) . (A9)After having performed here the Gaussian integrationover u and h , we are left with an integral over dχ , which3can be calculated by the stationary-phase method. De-noting the value of χ obeying the stationary-phase con-dition as χ = iσ , we express the free energy F in terms of σ and ξ : FL = − σ µ + λ ) + σ ( ξ −
1) (A10)+ 12 Z q ω (cid:8) d c ln (cid:0) κ q + σ q + ρω (cid:1) + ln (cid:2) κ q + ( σ + µξ ) q + ρω (cid:3) + ln (cid:2) κ q + ( σ + [2 µ + λ ] ξ ) q + ρω (cid:3)(cid:9) . Effects of the anharmonic coupling between in-planeand flexural modes lead to the following modifications ofEq. (A10). First, the bending rigidity κ gets renormal-ized, κ → κ q , as discussed in Sec. I. Second, there arisesa self-energy correction [37] σ to σ in the argumentsof logarithms. We will calculate σ in Appendix B. Asdiscussed in Sec. III A and in Appendix B, the total coef-ficient of the q term in the phonon propagator is exactlythe external tension σ , σ + σ = σ. (A11)This relation is a manifestation of a Ward identity thatis verified within the RG analysis (in one-loop order) inAppendix D.The stationary-point condition ∂F/∂σ = 0 yields σ µ + λ = ξ −
1+ 12 Z q ω (cid:20) d c q κ q +( σ + σ ) q + ρω + q ( σ + σ + µξ ) q + ρω + q [ σ + σ + (2 µ + λ ) ξ ] q + ρω (cid:21) . (A12)Here, we have neglected the term κ q which is smallcompared to µξ q and (2 µ + λ ) ξ q . The stretching pa-rameter ξ entering this equation is fixed by the external tension σ , which is given by a derivative of the free en-ergy with respect to the“projected area” of the sample, A = ξ L (see Eq. (15)). Substituting Eq. (A10) intoEq. (15), we get σ = σ + 12 Z q ω (cid:20) µq ( σ + σ + µξ ) q + ρω + (2 µ + λ ) q [ σ + σ + (2 µ + λ ) ξ ] q + ρω (cid:21) . (A13)Neglecting in the denominator the tension σ + σ = σ (which is assumed to be much smaller than the in-planeelastic moduli µ and λ ), we get σ = σ + δσ, (A14)where δσ = 12 Z q ω (cid:20) µq µξ q + ρω + (2 µ + λ ) q (2 µ + λ ) ξ q + ρω (cid:21) . (A15)We will show in Appendix B that the one-loop contri-butions σ and δσ arising in the analysis of the tensionon the basis of the FP Green function and of the ther-modynamic relation, Eq. (15), respectively, are identical, σ = δσ , so that Eqs. (A11) and (A14) are fully consis-tent. Appendix B: Anharmonic coupling between in-planeand out-of-plane phonons
In this Appendix, we provide a derivation of the self-energy correction σ in Eq. (A12), which is generated bythe anharmonic coupling between in-plane and out-of-plane modes. In combination with results of Appendix A,this allows us to derive equations (14) and (20) the maintext used there for the analysis of thermomechanicalproperties of the membrane.We begin with the action (A8). Integrating out thein-plane modes, we get an energy functional which onlydepends on h − fields: [37] E [ h ] = 12 Z q Ω ( κ q + ρω ) | h q ω | + 18 Z q Ω Z Q ω Z Q ′ ω ′ R αβγθ ( q , Ω)( Q α − q α ) Q β ( Q ′ γ + q γ ) Q ′ θ ( h ω Q h Ω − ω q − Q )( h ω ′ Q ′ h − Ω − ω ′ − q − Q ′ ) (B1)Here R αβγθ ( q , Ω) is the coupling tensor with the following nonzero components [in the basis of vectors n =4( − q y , q x ) / | q | , m = q / | q | ]: R nnnn = 4 µ ( µ + λ )2 µ + λ + ρ Ω λ (2 µ + λ )[ q (2 µ + λ ) ξ + ρ Ω ] , (B2) R mmmm = ρ Ω (2 µ + λ ) q (2 µ + λ ) ξ + ρ Ω , (B3) R mmnn = R nnmm = ρ Ω λq (2 µ + λ ) ξ + ρ Ω , (B4) R mnmn = 4 ρ Ω µq µξ + ρ Ω . (B5)The R αβγθ couplings characterize the quartic interactionof the out-of-plane modes. This interaction generatesself-energies the hh correlation functions. The couplingconstant R nnnn leads to a self-energy that scales as q ln q and, therefore, is responsible for the power-law renormal-ization of κ . On the other hand, the couplings R mmmm and R mnmn lead to self-energy corrections that scale as q . Specifically, the self-energy originating from the cou-pling R mmmm readsΣ mmmm Q , Ω = Z q ω ( Qq − q ) ( Qq ) q ρ (2 µ + λ ) ω (2 µ + λ ) q ξ + ρω × κ ( q − Q ) + ρ ( ω − Ω) (B6)(the factor 1 / h fields.) In the limit Ω → Q → , we getΣ mmmmQ → , Ω → = Q Z q ω ρ (2 µ + λ ) q ω (2 µ + λ ) ξ q + ρω κ q + ρω (B7)The integral is determined by the ultraviolet cut-off ofthe theory q uv . For q < q uv , one can neglect the term κ q in the denominator. This yieldsΣ mmmm ≃ Q σ mmmm , with σ mmmm ≈ Z q ω (2 µ + λ ) q (2 µ + λ ) ξ q + ρω . (B8)Proceeding in the same way, we find σ mnmn ≈ Z q ω µq µξ q + ρω . (B9)Combining these contributions, we get the following re-sult for the total one-loop coefficient σ = σ mmmm + σ mnmn of the Q self-energy: σ ≈ Z q ω q (cid:20) µ + λ (2 µ + λ ) ξ q + ρω + µµξ q + ρω . (cid:21) (B10)Comparing this equation with Eq. (A15), we satisfy our-selves that σ = δσ, (B11)as was stated in the end of Appendix A. In combinationwith Eq. (A14), this yields Eq. (A11).We have thus explicitly demonstrated that Eq. (A12)can be written in terms of the applied tension σ , σµ + λ = ξ −
1+ 12 Z q ω d c q κ q + σq + ρω + a, (B12)where a = 12( µ + λ ) Z q ω q (cid:20) µ + λ ( σ + µξ ) q + ρω (B13)+ 3 µ + 2 λ [ σ + (2 µ + λ ) ξ ] q + ρω (cid:21) is an ultraviolet correction that can be fully absorbed inthe renormalization of ξ. Equation (B12) is the gener-alized Hooke’s law. We emphasize once more that thedenominator in the r.h.s. of Eq. (B12) contains only theexternal tension σ and is not sensitive to the ultravioletcutoff of the theory.Since a depends on T, it leads to a correction to α T . One can show that this correction is of the order of δα T ∼ κ (cid:18) Tg κ (cid:19) , and is thus small in comparison with Eq. (38) under thecondition T < g κ p /η + ln(1 /g ) (for graphene T <
T < T uv , with T uv given by Eq. (28), which is the temperature range of ourinterest in this paper. Neglecting a in Eq. (B12), we getEq. (20) of the main text.5Using the identity (A11), we can also cast the free en-ergy of renormalized out-of-plane modes in Eq. (A10) inthe form FL = − σ B + σ ( ξ −
1) (B14)+ d c X q ω ln (cid:0) κ q q + σq + ρω (cid:1) , which is Eq. (14) of the main text. Appendix C: Quantum renormalization group.
In this Appendix, we derive the quantum RG equationsthat are presented in Sec. III B of the main text. Analternative derivation is presented in Appendix D.The quantum renormalization grouo operates in theregion of momenta q T < q < q uv . For momenta q be-low q uv ∼ p µ/ κ , the flexural phonons are softer thanthe in-plane modes: ω q < ω ⊥ , k q . Here ω ⊥ q = p µ/ρq and ω k q = p (2 µ + λ ) /ρq . Since in the considered re-gion of momenta ~ ω q ≫ T , the flexural phonons arefrozen out, so that the relevant RG equations are of zero-temperature character [62].Here, we derive RG equations by using the energy func-tional Eq. (B1) where in-plane modes have been inte-grated out.We have demonstrated above that the terms scalingas q in the flexural phonons self-energy cancel (for zeroexternal tension σ = 0). After this cancellation is takeninto account, all remaining effects related to retardationturn out to be small for q < q uv and can be safely ne-glected. Hence, it is sufficient to keep the only componentof the interaction tensor R αβγθ ,R nnnn ≈ µ ( µ + λ )2 µ + λ = Y. (C1)To proceed, we use the approach analogous to one de-veloped in Ref. 39 for high-temperature case. To findthe renormalization of elastic coefficients within this ap-proach, we have to calculate the polarization operatorand the self-energy of h − fields. The bare Green functionfor h -field reads G ω, k = 1 κ k + ρω . (C2)The polarization operator is given by the following equa-tion Π Ω , q = d c Z dωd k (2 π ) k ⊥ G ω, k , G − ω, q − k , (C3)where k ⊥ = [ k × q ] /q . Equation (C3) is a quantum coun-terpart of Eq. (37) of Ref. 39 derived there for the high-temperature classical regime. Performing the integration over dω , we getΠ Ω , q = d c ρ Z d k (2 π ) k ⊥ ( ω k − q + ω k ) ω k − q ω k [( ω k − q + ω k ) + Ω ] (C4)Carrying out the remaining momentum integration wefind, in the limit Ω = 0 and q → q = Π Ω=0 , q → ≈ d c πρ / κ / ln (cid:18) q uv q (cid:19) . (C5)Next, we use Eqs. (33) and (34) of Ref. 39 to find screen-ing of coupling constants Y and µ : Y q = Y Y Π q / ≈ Y − Y Π q / , (C6) µ q = µ µ Π q ≈ µ − µ Π q . (C7)We notice that for D = 2 the interaction be-tween h − fields in Eq. (B1) depends on coupling Y = R nnnn (Ω →
0) only, while µ drops out from Eq. (B1).Therefore, in order to obtain Eq. (C7), one should firstconsider D = 2 , and then take the limit D → . Substituting Eq. (C5) into Eqs. (C6) and (C7), we findRG equations for in-plane elastic moduli: dYd
Λ = − d c Y πρ / κ / , (C8) d ( µ + λ ) d Λ = − d c ( µ + λ ) πρ / κ / , (C9)where Λ = ln ( q uv /q ) . These equations are equivalentto Eqs. (9) and (10) of Ref. 62. Renormalization ofself-energy of h − field is given by an equation similar toEq. (43) of Ref. 39:Σ ω, k = Z d Ω d q (2 π ) k ⊥ Y q G ω − Ω , k − q . (C10)Integrating over d Ω , taking the limit ω → , and neglect-ing the dependence of Y on q, we getΣ ω → , k = Y √ κ ρ Z d q (2 π ) k ⊥ | k − q | . (C11)A straightforward analysis of this integral shows that Σscales as k ln k, which implies a renormalization of κ , d κ d Λ = 3 Y π √ κ ρ . (C12)The Eq. (C12) coincides up to the sign with Eq. (11)of Ref. 62. From Eqs. (C6),(C7) and (C12), one easilyobtains Eqs. (29) and (30) of the main text, with g givenby Eq. (12).6 Appendix D: Background-field renormalization
In this Appendix, we perform a derivation of quantum RG equations (Sec. III B of the main text) alternative to thatpresented in Appendix C. For this purpose, we evaluate the self-energies of the propagators of in- and out-of-planephonons within one-loop approximation by using the approach of Ref. [62].We start from the Lagrangian given by Eq. (A2). Using definition (15), we obtain formally exact relation betweenthe external tension and the global stretching factor: σB = ξ − Z q ,ω q d c G ω,q + 2 µ + λ B F ( t ) ω,q + 3 µ + 2 λ B F ( l ) ω,q ! + 12 ξB ( µδ γβ δ ηα + λδ γα δ ηβ ) (cid:10) ∂ α u γ (cid:0) ∂ η u ∂ β u + ∂ η h ∂ β h (cid:1)(cid:11) . (D1)Here F ( t,l ) ω,q and G ω,q are exact (with respect to the full Lagrangian (A2)) propagators of in- and out-of plane phonons,respectively: h u α ( q , iω ) u β ( − q , − iω ) i = F ( l ) ω,q q α q β q + F ( t ) ω,q (cid:18) δ αβ − q α q β q (cid:19) , h h ( q , iω ) h ( − q , − iω ) i = d c G ω,q . (D2)The exact propagators can be cast in the following form[ F ( l ) ω,q ] − = ρω + [(2 µ + λ ) ξ + ( µ + λ )( ξ − q + κ q − Σ ( l ) ω,q , [ F ( t ) ω,q ] − = ρω + [ µξ + ( µ + λ )( ξ − q + κ q − Σ ( t ) ω,q , [ G ω,q ] − = ρω + ( µ + λ )( ξ − q + κ q − Σ ω,q , (D3)where the self-energies Σ ( l,t ) ω,q and Σ ω,q take into account interaction of in- and out-of-plane modes encoded in theLagrangian (A5). The expressions (D3) in the absence of self-energies corresponds to the Gaussian part of theLagrangian (A2). We emphasize the appearance of linear in q term in the propagator of the out-of-plane phonondue to the linear in K term in the Lagrangian (A2). As we will see below it will be compensated by the linear in q term from the self-energy Σ ω,q . To avoid confusion, we note that the definition of the self-energies used in AppendixesA, B, and C is different compared to the definition which we use here. Of course, this does not change the physicalpropagators and, in particular, the cancelation of ∝ q contributions to the inverse propagator of out-of-plane phononsin the absence of the external stress ( σ = 0). This statement can be written as σ + σ = 0 (as was done in AppendixesA, B, C), or, equivalently, as B ( ξ − − lim q → Σ ω =0 ,q /q = 0 within background-filed renormalization approach usedin this Appendix.In order to find the corresponding self-energies Σ ( l,t ) ω,q and Σ ω,q , we use the background field method. We split thefields u and h on slow u ′ , h ′ and fast ˜ u , ˜ h components in the momentum and frequency spaces, u = u ′ + ˜ u and h = h ′ + ˜ h . We denote the corresponding momentum scale which separates fast and slow modes as q sfΛ . Then theinteraction terms in the Lagrangian (A5) generates the following interaction terms between slow and fast components.For a sake of simplicity, we consider the case of d c = 1 and restore arbitrary dimensionality in the final results for the7self-energies only. Then, limiting ourselves to the first and second orders in slow components, we find S (1) , u ′ , ˜ h = Z β dτ Z d x (cid:16) µξ∂ α u ′ β + λξ δ αβ ∂ η u ′ η (cid:17) ∂ α ˜ h∂ β ˜ h,S (1) , , h ′ , ˜ u , ˜ h = Z β dτ Z d x h µξ (cid:16) ∂ α ˜ u β + ∂ β ˜ u α (cid:17) + λξδ αβ ∂ η ˜ u η i ∂ α h ′ ∂ β ˜ h,S (1) , u ′ , ˜ u = Z β dτ Z d x h(cid:16) µξ∂ α u ′ β + λξ δ αβ ∂ η u ′ η (cid:17) ∂ α ˜ u γ ∂ β ˜ u γ + (cid:16) µξ (cid:16) ∂ α ˜ u β + ∂ β ˜ u α (cid:17) + λξδ αβ ∂ η ˜ u η (cid:17)i ∂ α u ′ γ ∂ β ˜ u γ ,S (2) , u ′ , ˜ u = 12 (cid:16) µδ αθ δ βη + λ δ αη δ βθ (cid:17) Z β dτ Z d x (cid:16) ∂ α u ′ γ ∂ η u ′ γ ∂ θ ˜ u ζ ∂ β ˜ u ζ + ∂ α u ′ γ ∂ η ˜ u γ ∂ θ ˜ u ζ ∂ β u ′ ζ + ∂ α u ′ γ ∂ η ˜ u γ ∂ θ u ′ ζ ∂ β ˜ u ζ (cid:17) ,S (2) , h ′ , ˜ h = 14 (cid:16) µδ αθ δ βη + λ δ αη δ βθ (cid:17) Z β dτ Z d x (cid:16) ∂ α h ′ ∂ η h ′ ∂ θ ˜ h∂ β ˜ h + (cid:0) ∂ α h ′ ∂ η ˜ h + ∂ α ˜ h∂ η h ′ (cid:1)(cid:0) ∂ θ ˜ h∂ β h ′ + ∂ θ h ′ ∂ β ˜ h (cid:1)(cid:17) ,S (2) , u ′ , ˜ h = 12 (cid:16) µδ αθ δ βη + λ δ αη δ βθ (cid:17) Z β dτ Z d x ∂ α u ′ γ ∂ η u ′ γ ∂ θ ˜ h∂ β ˜ h,S (2) , h ′ , ˜ u = 12 (cid:16) µδ αθ δ βη + λ δ αη δ βθ (cid:17) Z β dτ Z d x ∂ α h ′ ∂ η h ′ ∂ θ ˜ u γ ∂ β ˜ u γ ,S (1 , , , h ′ , u ′ , ˜ h, ˜ u = 12 (cid:16) µδ αθ δ βη + λ δ αη δ βθ (cid:17) Z β dτ Z d x (cid:0) ∂ α u ′ γ ∂ η ˜ u γ + ∂ α ˜ u γ ∂ η u ′ γ (cid:1) ( ∂ θ ˜ h∂ β h ′ + ∂ θ h ′ ∂ β ˜ h (cid:1) . (D4)After integration over fast variables we find the correction to the Gaussian part of the action for the in-plane ( δS (2) u )and out-of-plane ( δS (2) h ) slow modes: δS (2) u = h S (2) , u ′ , ˜ h i + h S (2) , u ′ , ˜ u i − (cid:28)h S (1) , u ′ , ˜ u i (cid:29) − (cid:28)h S (1) , u ′ , ˜ h i (cid:29) , (D5) δS (2) h = h S (2) , h ′ , ˜ h i + h S (2) , h ′ , ˜ u i − (cid:28)h S (1) , , h ′ , ˜ u , ˜ h i (cid:29) . (D6)Here the average h . . . i is with respect to the Gaussian part of the full Lagrangian (A2). The self-energies can befound from the following expressions δS (2) u = − Z q ,ω u ′ α ( q , iω ) u ′ β ( − q , − iω ) (cid:20) Σ ( l ) ω,q q α q β q + Σ ( t ) ω,q (cid:18) δ αβ − q α q β q (cid:19)(cid:21) ,δS (2) h = − Z q ,ω h ′ ( q , iω ) h ′ ( − k , − iω )Σ ω,q . (D7)Evaluation of averages in δS (2) h yieldsΣ ω,q = − q Z k , Ω " (2 µ + λ ) k G (0)Ω ,k + λ + µ k (cid:16) F ( l ) , (0)Ω ,k + F ( t ) , (0)Ω ,k (cid:17) − Z k , Ω ξ G (0) ω +Ω , q + k nh (2 µ + λ ) k ( q · k ) +4 µ (2 µ + λ )( q · k ) + 2 λ (2 µ + λ )( q · k ) q k + 4 µ ( q · k ) k − + 4 µλ ( q · k ) q + λ k q i F ( l ) , (0)Ω ,k + µ [ k × q ] k (cid:0) k + 2( k · q ) (cid:1) F ( t ) , (0)Ω ,k o . (D8)Here F ( l,t ) , (0)Ω ,k and G (0)Ω ,k denote the propagators of in- and out-of-plane modes within the Gaussian approximation tothe full Lagrangian (A2). In Eq. (D8) the terms linear in the propagators corresponds to the contributions from h S (2) , h ′ , ˜ h i and h S (2) , h ′ , ˜ u i whereas the terms proportional to the product of propagators for the in-plane and out-of-planemodes corresponds to the last contribution in the right hand side of Eq. (D6). We note that the first and last terms inthe right hand side of Eq. (D8) corresponds to the self-energy contribution due to effective interaction tensor R αβγθ .The second term in the right hand side of Eq. (D8) appears due to the interaction of two flexural phonons with two8 TABLE I: The linear in q contributions to Σ ( l,t ) ω =0 ,q from different terms in Eq. (D5). − Σ ( t ) ω =0 ,q /q − (Σ ( l ) ω =0 ,q − Σ ( t ) ω =0 ,q ) /q h S (2) , u ′ , ˜ h i µ + λ R k , Ω k G (0)Ω ,k h S (2) , u ′ , ˜ u i R k , Ω k (cid:16) µ +5 λ F ( l ) , (0)Ω ,k + µ +7 λ F ( t ) , (0)Ω ,k (cid:17) µ + λ R k , Ω k (cid:16) F ( l ) , (0)Ω ,k − F ( t ) , (0)Ω ,k (cid:17) − (cid:28)h S (1) , u ′ , ˜ u i (cid:29) − R k , Ω k ξ F ( l ) , (0)Ω ,k (cid:16) (3 µ + λ ) F ( l ) , (0)Ω ,k − R k , Ω k ξ (cid:16) (3 µ +2 λ ) [ F ( l ) , (0)Ω ,k ] + (2 µ + λ ) [ F ( t ) , (0)Ω ,k ] + µ +10 µλ +3 λ F ( t ) , (0)Ω ,k (cid:17) − ( µ + λ ) F ( l ) , (0)Ω ,k F ( t ) , (0)Ω ,k (cid:17) − (cid:28)h S (1) , u ′ , ˜ h i (cid:29) − µ ξ R k , Ω k [ G (0)Ω ,k ] − ( µ + λ ) ξ R k , Ω k [ G (0)Ω ,k ] in-plane phonons. This vertex is not included in the interaction tensor R αβγθ . However, as we shall demonstratebelow, this interaction is taken into account in the approach of Appendix A.In the limit q → ω = 0 we find from Eq. (D8)Σ ω =0 ,q = − q Z k ,ω " (2 µ + λ ) k G (0)Ω ,k + λ + µ k (cid:16) F ( l ) , (0)Ω ,k + F ( t ) , (0)Ω ,k (cid:17) + q ξ Z k , Ω k G (0)Ω ,k h (2 µ + λ ) F ( l ) , (0)Ω ,k + µ F ( t ) , (0)Ω ,k i . (D9)It is convenient to regroup various terms in Eq. (D9) in the following way:Σ ω =0 ,q = − q µ + λ Z k ,ω k " G (0)Ω ,k + F ( l ) , (0)Ω ,k + F ( t ) , (0)Ω ,k + q Z k , Ω k G (0)Ω ,k h (2 µ + λ ) ξ k F ( l ) , (0)Ω ,k + µ ξ k F ( t ) , (0)Ω ,k − (3 µ + λ ) i . (D10)As one can check, Eq. (D10) can be written as Σ ω =0 ,q = [ σ + σ − B ( ξ − q . Using the precise form of theGaussian propagators the result (D10) can be equivalently rewritten as followsΣ ω =0 ,q = − q Z k , Ω k ( d c ( µ + λ ) G (0)Ω ,k + (3 µ + 2 λ ) F ( l ) , (0)Ω ,k + (2 µ + λ ) F ( t ) , (0)Ω ,k ) . (D11)Here we restore arbitrary value of d c . Comparing this result with the expression (D1) evaluated within the Gaussiantheory, we conclude that within one-loop approximation the following identity holds σ = B ( ξ − − lim q → Σ ω =0 ,q /q . (D12)Although at present we cannot prove this relation beyond the one-loop approximation, we believe that it should besatisfied in general (see discussion in the main text).Expansion of the self-energy (D8) to the second in q determines the one-loop renormalization of the bendingrigidity: κ ′ = κ − ∂ ∂q Σ ω =0 ,q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) q =0 . (D13)As one can check by inspection of various terms in Eq. (D8), ∂ Σ ω =0 ,q /∂q , the logarithmically divergent contributionsappear only for external momentum scale q T < q < q uv ∼ p min { λ, µ } / κ . Simplifying Eq. (D8) in this regime, weobtain Eq. (C11) where the integration over momentum is limited to k > q sfΛ whereas q ≪ q sfΛ . Performing integrationover momentum, we find the following RG equation: d κ d Λ = 3 d c πρ / κ / µ ( µ + λ )2 µ + λ , (D14)where Λ = ln q uv /q since the minimal value of q sfΛ is given by the external momentum q . This equation coincides withEq. (30) in the main text.9The self-energies Σ ( t ) ω,q and Σ ( l ) ω,q determine renormalization of the Lame coefficients: µ ′ = µ − lim q → h Σ ( t ) ω =0 ,q − Σ ω =0 ,q i. ( ξ q ) , ( λ + µ ) ′ = λ + µ − lim q → h Σ ( l ) ω =0 ,q − Σ ω =0 ,q i. ( ξ q ) . (D15)We do not present the full expressions for the self-energies Σ ( l,t ) ω,q since they are too cumbersome. The linear in k contributions to Σ ( l,t ) ω =0 ,q from different terms in Eq. (D5) are summarized in Table I. Using Eq. (D15) and Table I,we find µ ′ = µ − Z k , Ω k ( d c µ (cid:2) G (0)Ω ,k (cid:3) + (3 µ + λ ) (cid:2) F ( l ) , (0)Ω ,k ] + 2 µ ( λ + 2 µ ) F ( l ) , (0)Ω ,k F ( t ) , (0)Ω ,k ) , (D16)( λ + µ ) ′ = λ + µ − Z k , Ω k ( d c ( µ + λ ) (cid:2) G (0)Ω ,k (cid:3) + (3 µ + 2 λ ) (cid:2) F ( l ) , (0)Ω ,k ] + (2 µ + λ ) (cid:2) F ( t ) , (0)Ω ,k ] ) . (D17)Here we restore arbitrary value of d c .Assuming that the infrared moment scale (which sep-arates the slow and fast modes in the moment space) liesin the range q T < q < q uv ∼ p min { λ, µ } / κ we find thatonly the terms proportional to d c provide logarithmicallydivergent contributions in Eqs. (D16) and (D17). Hence,we find dµd Λ = − d c πρ / κ / µ ,dλd Λ = − d c πρ / κ / ( µ + 4 λµ + 2 λ ) . (D18)From Eqs. (D14) and (D18) we obtain the renormal-ization group equations (29) and (30) of the main text.We see that q uv ∼ p min { λ, µ } / κ is a natural ultravioletcut-off for the renormalization group equations (29) and(30).As we mentioned in the main text, the momentum q uv ∼ p min { λ, µ } / κ is on the order of the inverse latticeconstant a − for graphene. However, one can imagine a generic membrane, where q uv ≪ /a. Let us briefly dis-cuss what happens for q uv < q < Q uv ∼ /a. Withinthis interval, there is no difference in the spectrum of in-plane and out-of-plane phonons. Then all terms in theright hand side of Eqs. (D16) and (D17) provide loga-rithmic contributions. Then we find the following renor-malization group equations for the Lame coefficients inthe range q uv < q < Q uv : dµd ˜Λ = − ( λ + 4 µ ) + ( d c − µ πρ / κ / ,dλd ˜Λ = − d c ( λ + µ ) + (3 λ + 2 µ ) + (7 − d c ) µ πρ / κ / . (D19)where ˜Λ = ln Q uv /q . We note that in the range q uv < q Encyclopedia of Condensed Matter Physics (Elsevier,Amsterdam, 2005), ed. by G. F. Bassani, G. L. Liedl, andP. Wyder, p. 77.[7] W. Miller, C. W. Smith, D. S. MacKenzie, K. E. Evans,Journal of Materials Science , 241402(R) (2001)[9] D. Tomanek, J. Phys.: Condens. Matter R413 (2005).[10] Y.-K. Kwon, S. Berber, and D. Tomanek, Phys. Rev.Lett. , 015901 (2004).[11] S. Brown, J. Cao, J. L. Musfeldt, N. Dragoe, F. Cim-poesu, S. Ito, H. Takagi, and R. J. Cross, Phys. Rev. B , 125446 (2006).[12] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y.Zhang, S. V. Dubonos, I. V. Grigorieva and A.A. Firsov,Science , 666 (2004).[13] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,M. I. Katsnelson, I.V. Grigorieva, S.V. Dubonos and A.A. Firsov, Nature , 197 (2005).[14] Y. Zhang, Y.-W. Tan, H. L. Stormer and P. Kim, Nature , 201 (2005).[15] A. K. Geim and K. S. Novoselov, Nature Materials , 183 (2007).[16] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009).[17] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi,Rev. Mod. Phys. , 407 (2011).[18] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, andA. H. Castro Neto, Rev. Mod. Phys. , 1067 (2012).[19] M. I. Katsnelson, Graphene: Carbon in Two DimensionsHardcover , Cambridge University Press (2012).[20] E. L. Wolf, Graphene: A New Paradigm in CondensedMatter and Device Physics , Oxford University Press(2014).[21] L. E. F. Foa Torres, S. Roche, J.-C. Charlier, Introduc-tion to Graphene-Based Nanomaterials From ElectronicStructure to Quantum Transport , Cambridge UniversityPress (2014).[22] M. K. Blees, A.W. Barnard, P. A. Rose, S. P. Roberts,K. L. McGill, P. Y. Huang, A. R. Ruyack, J. W. Kevek,B. Kobrin, D. A. Muller, and P. L. McEuen, Nature ,204 (2015).[23] G. Lopez-Polin, C. Gomez-Navarro, V. Parente, F.Guinea, M. I. Katsnelson, F. Perez-Murano, and J.Gomez-Herrero, Nature Physics , 26 (2015); G. Lopez-Polin, M. Jaafar, F. Guinea, R. Roldan, C. Gomez-Navarro, and J. Gomez-Herrero, arXiv:1504.05521.[24] R. J. T. Nicholl, H. J. Conley, N. V. Lavrik, I. Vlas-siouk, Y. S. Puzyrev, V. P. Sreenivas, S. T. Pantelides,and K. I. Bolotin, Nature Communications 6:8789 doi:10.1038/ncomms9789 (2015).[25] I. M. Lifshitz, Sov. Phys. JETP , 472 (1952).[26] N. Mounet and N. Marzari, Phys. Rev. B , 205214(2005).[27] P. L. de Andres, F. Guinea, and M. I. Katsnelson, Phys.Rev. B , 144103 (2012).[28] M. I. Katsnelson and A. Fasolino, Acc. Chem. Res. ,97 (2013).[29] E. G. Steward, B. P. Cook, and E. A. Kellert, Nature , 046808 (2009).[31] A. L. C. da Silva, Ladir Cˆandido, J. N. Teixeira Ra-belo, G.-Q. Hai, and F. M. Peeters, Europhys. Lett., ,56004 (2014).[32] K. H. Michel, S. Costamagna, and F. M. Peeters, Phys.Rev. B , 134302 (2015).[33] W. Bao, F. Miao, Z. Chen, H. Zhang, W. Jang, C. Dames,and C. N. Lau, Nature Nanotech. 4, 562.[34] D. Yoon, Y.-W. Son, and H. Cheong, Nano Lett. 11, 3227.[35] F. Boerrnert, A. Barreiro, D. Wolf, M. I. Katsnelson, B.Buechner, L. M. K. Vandersypen, and M. H. Ruemmeli,Nano Lett. , 4455 (2012).[36] L. D. Landau and E. M. Lifshitz, Statistical Physics, PartI (Pergamon Press, Oxford, 1980).[37] B. Amorim, R. Roldan, E. Cappelluti, A. Fasolino, F.Guinea, and M. I. Katsnelson, Phys. Rev. B , 224307(2014).[38] D. Nelson, T. Piran, and S. Weinberg (Eds.) StatisticalMechanics of Membranes and Surfaces (World Scientific,Singapore, 1989).[39] I. V. Gornyi, V. Yu. Kachorovskii, and A. D. Mirlin,Phys. Rev. B , 155428 (2015).[40] D. R. Nelson and L. Peliti, J. Phys. (Paris) , 1085 (1987).[41] Y. Kantor and D. R. Nelson, Phys. Rev. Lett. , 2774(1987); Phys. Rev. A , 4020 (1987);[42] M. Paczuski, M. Kardar, and D. R. Nelson, Phys. Rev.Lett. , 2638 (1988).[43] F. David and E. Guitter, Europhys. Lett. , 709 (1988).[44] E. Guitter, F. David, S. Leibler, and L. Peliti, Phys.Rev. Lett. , 2949 (1988).[45] J. A. Aronovitz and T. C. Lubensky, Phys. Rev. Lett. , 2634 (1988).[46] E. Guitter, F. David, S. Leibler, and L. Peliti, J.Phys.France 609 (1989).[48] M. Paczuski and M. Kardar, Phys. Rev. A , 6086(1989).[49] L. Radzihovsky and D. R. Nelson, Phys. Rev. A , 3525(1991).[50] D. R. Nelson and L. Radzihovsky, Europhys. Lett. ,79 (1991).[51] G. Gompper and D. M. Kroll, Europhys. Lett. , 783(1991).[52] L. Radzihovsky and P. Le Doussal, J.Phys. I France 599 (1992).[53] P. Le Doussal and L. Radzihovsky, Phys. Rev. Lett ,1209 (1992).[54] D. C. Morse, T. C. Lubensky, and G. S. Grest, Phys.Rev. A , R2151 (1992).[55] P. Le Doussal and L. Radzihovsky, Phys. Rev. B , 3548(1993).[56] M. J. Bowick, S. M. Catterall, M. Falcioni, G. Thorleif-sson, and K. N. Anagnostopoulos, J. Phys. I France ,1321 (1996).[57] I.V. Gornyi, V. Yu. Kachorovskii, and A. D. Mirlin, 2DMaterials , 011003 (2016).[58] J.-P. Kownacki, and D. Mouhanna, Phys. Rev. E ,040101(R) (2009).[59] D. Gazit, Phys. Rev. E , 041117 (2009).[60] F. L. Braghin and N. Hasselmann , Phys. Rev. B ,035407 (2010).[61] V. V. Lebedev and E. I. Kats, Phys. Rev. B , 045416(2012).[62] E. I. Kats and V. V. Lebedev, Phys. Rev. B , 125433(2014).[63] E. I. Kats and V. V. Lebedev, Phys. Rev. B , 176301(2014).[64] E. I. Kats and V. V. Lebedev, Phys. Rev. B ,079904(E) (2016).[65] J. H. Los, M. I. Katsnelson, O. V. Yazyev, K. V. Za-kharchenko, and A. Fasolino, Phys. Rev. B , 121405(2009).[66] E. I. Kats and V. V. Lebedev, Phys. Rev. E , 032415(2015).[67] B. Amorim, R. Roldan, E. Cappelluti, F. Guinea, A.Fasolino, and M. I. Katsnelson, Phys. Rev. B , 176302(2014)[68] O. Coquand, D. Mouhanna, unpublished, cond-mat arxiv1607.03335[69] P. Souvatzis, S. Arapan, O. Eriksson, and M. I. Katsnel-son, Europhys. Lett. , 66006 (2011).[70] Analyzing integral Eq. (46), one can find that its low σ asymplotics has form of the upper line in Eq. (48) pro- vided that C − η ) = (cid:20) A ( η )2 − η (cid:21) − η (cid:18) B Y (cid:19) η , where A ( η ) = Z ∞ ηdyy − η (1 + y − η ) = 2 − η (cid:18) − η − η (cid:19) Γ (cid:18) − η − η (cid:19) . We notice, that for invariant subspace of elastic moduli, Y = 8 B/ , the coefficient C is temperature-independent.For d c → ∞ ( η → C = 1 / . [71] J. H. Los, A. Fasolino, and M. I. Katsnelson, Phys. Rev.Lett. , 015901 (2016).[72] A. K. Geim and I. V. Grigorieva, Nature , 419 (2013).[73] K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H.Castro Neto, Science , 461 (2016). [74] C. R.Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C.Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gor-bachev, A. V. Kretinin, J. Park, L. A. Ponomarenko,M. I. Katsnelson, Yu. N. Gornostyrev, K. Watanabe, T.Taniguchi, C. Casiraghi, H.-J. Gao, A. K. Geim, and K.S. Novoselov, Nature Phys. , 451 (2014).[75] M. M. van Wijk, A. Schuring, M. I. Katsnelson, and A.Fasolino, Phys. Rev. Lett. , 135504 (2014).[76] C. R. Woods, F. Withers, M. J. Zhu, Y. Cao, G.Yu,A. Kozikov, M. Ben Shalom, S. V. Morozov, M. M.van Wijk, A. Fasolino, M. I. Katsnelson, K.Watanabe,T.Taniguchi, A. K. Geim, A. Mishchenko, and K. S.Novoselov, Nature Communications , 10800 (2016)[77] E. Khestanova, F. Guinea, L. Fumagalli, A. K. Geim,and I. V. Grigorieva, Nature Commun.7