The effect of anomalous elasticity on the bubbles in van der Waals heterostructures
TThe effect of anomalous elasticity on the bubbles in van der Waals heterostructures
A. A. Lyublinskaya, S. S. Babkin, and I. S. Burmistrov
2, 3 Moscow Institute of Physics and Technology, 141700, Dolgoprudnyi, Moscow Region, Russia L. D. Landau Institute for Theoretical Physics, Semenova 1-a, 142432, Chernogolovka, Russia Laboratory for Condensed Matter Physics, National Research University Higher School of Economics, 101000 Moscow, Russia
It is shown that the anomalous elasticity of membranes affects the profile and thermodynamics of abubble in van der Waals heterostructures. Our theory generalizes the non-linear plate theory as wellas membrane theory of the pressurised blister test to incorporate the power-law scale dependence ofthe bending rigidity and Young’s modulus of a two-dimensional crystalline membrane. This scaledependence caused by long-ranged interaction of relevant thermal fluctuations (flexural phonons),is responsible for the anomalous Hooke’s law observed recently in graphene. It is shown that thisanomalous elasticity affects dependence of the maximal height of the bubble on its radius andtemperature. We identify the characteristic temperature above which the anomalous elasticity isimportant. It is suggested that for graphene-based van der Waals heterostructures the predictedanomalous regime is experimentally accessible at the room temperature.
Introduction. — Mechanical properties of two-dimensional (2D) materials, especially, of van der Waalsheterostructures, have recently attracted a lot of interestin view of their potential applications [1]. The simplestexample of van der Waals heterostructure is two mono-layers, e.g. graphene, hexagonal boron nitride (hBN),MoS , assembled together. Strong adhesion betweenmonolayers [2] results in atomically clean interfaces inwhich all contaminating substances are combined intobubbles [3]. Recently, these bubbles inside van der Waalsheterostructures have been studied experimentally [4].Similar bubbles are formed between an atomic monolayerand a substrate, e.g. SiO [4, 5]. There are many sug-gestions of practical usage of the bubbles inside van derWaals heterostructures, for example, the graphene liquidcell microscopy [6], controlled room-temperature photo-luminescence emitters [7], etc.The mechanics of monolayers due to these bubbles isconsidered to be analogous to the one of the pressur-ized blister test which has recently become the routinemethod to measure simultaneously the Young’s modu-lus and adhesion energy of a monolayer on a substrate[8–11]. Usually, the pressurized blister test is describedeither by nonlinear plate model or by membrane the-ory (see e.g. [12, 13]). These standard elastic theoriesof deformed plates ignore the fact that elastic proper-ties of an atomic monolayer are those of 2D crystallinemembranes [14–20]. The most striking feature of me-chanics of membranes is anomalous elasticity which re-sults in non-linear (so-called, anomalous) Hooke’s law forsmall applied stress (see Refs. [21, 22] for a review). Re-cently, this anomalous Hooke’s law has been measured ingraphene [23, 24]. However, until present, see e.g. Refs.[25–28], the anomalous elasticity of membranes is com-pletely ignored in description of mechanical properties ofmonolayers in the presence of the bubbles.The present paper generalizes the classical theory ofthe pressurised blister test to incorporate the anomalouselasticity of a membrane. Our approach explicitly takes FIG. 1. Sketch of a spherical bubble between a membraneand a substrate. into account the power-law renormalization of the bend-ing rigidity and Young’s modulus. It is shown that abovea certain temperature the dependence of the bendingrigidity and Young’s modulus of a membrane on the ra-dius of the bubble results in the non-analytic dependenceof its maximal height on the radius and temperature.
Model. — The model for the description of the profileof a bubble between the membrane and the substrate (seeFig. 1) is well established [4, 12]. It can be formulatedin terms of the (free) energy which is the sum of thefollowing four terms: E = E bend + E el + E b + E vdW . (1)Here the first contribution describes the energy cost re-lated with bending of the membrane: E bend = κ (cid:90) d r (cid:2) (∆ h ) + (∆ u ) (cid:105) , (2)where κ denotes the bare bending rigidity of the mem-brane. Here u = { u x , u y } and h are the in-plane andout-of-plane displacements of the membrane (see Fig. 1).The second term is the standard elastic energy [29]: E el = (cid:90) d r (cid:0) µ u αβ u βα + λ u αα u ββ / (cid:1) , (3)where µ and λ stand for the Lam´e coefficients and u αβ = (cid:0) ∂ β u α + ∂ α u β + ∂ α h∂ β h + ∂ α u ∂ β u (cid:1) / E b describes the con-tribution of the bubble substance. Under assumption a r X i v : . [ c ond - m a t . m e s - h a ll ] S e p of the constant pressure P inside the bubble one has E b = − P V , where the bubble volume can be approxi-mated as V = (cid:82) d r h ( r ). The last term in Eq. (1) de-scribes the van der Waals interaction between the mem-brane and the substrate in the presence of the bub-ble. It can be approximated as E vdW = πγR . Here R denotes the radius of the bubble (see Fig. 1) and γ = γ ms − γ mb − γ sb where γ ms , γ mb , and γ sb are theadhesion energies between the membrane and substrate,between the membrane and the substance inside the bub-ble, and between the substrate and the substance, respec-tively. The form (1) of the total energy is well justifiedif the maximal height of the bubble, H = h (0), is smallcompared to the radius, H (cid:28) R . Throughout the paperwe shall assume that this condition is fulfilled. The standard approach. — In order to compute the pro-file of the spherical bubble, initially, one needs to solvethe Euler-Lagrange equations for u ( r ) and h ( r ) with theproper boundary conditions and compute the energy E asa function of R and H . Usually, instead of the solutionof the Euler-Lagrange equations the approximate solu-tions either within the non-linear plate theory or withinthe membrane theory are used, see e.g. Ref. [12]. Fi-nally, one have to minimize E with respect to the both R and H . The minimization procedure allows one to findthe maximal height H and the pressure P as a functionof the bubble radius R . Comparison of the linear in u term and the term quadratic in h in the strain tensor u αβ leads to the following relation for the maximal hor-izontal deformation: u max ∼ H /R . In the consideredregime, H/R (cid:28)
1, the horizontal displacement is alsosmall, u max (cid:28) H . This implies that one can neglect theterm ∂ α u ∂ β u in u αβ . Also this allows one to omit theterm (∆ u ) in the bending energy such that it reads: E bend = κ (cid:90) d r (∆ h ) . (4)In the absence of ∂ α u ∂ β u in u αβ the elastic energybecomes quadratic in u . This implies that the Euler-Lagrange equation for u ( r ) become linear and can besolved for arbitrary configuration of h ( r ) (even not nec-essarily obeying the Euler-Lagrange equation). In otherwords, a horizontal deformation is adjusted to any verti-cal displacement. Therefore, E el is given as [14]: E el = Y (cid:90) d r (cid:20) K αα − ∂ α (cid:90) d r (cid:48) G ( r , r (cid:48) ) ∂ β K αβ ( r (cid:48) ) (cid:21) , (5)where K αβ = ∂ α h∂ β h and Y = µ ( µ + λ )2 µ + λ is the Young’smodulus. The function G ( r , r (cid:48) ) is the Green’s functionof the Laplace operator on the disk r (cid:54) R .Using Eqs. (4) and (5), we can estimate the bend-ing and elastic energies as E bend ∼ κ H /R and E el ∼ Y H /R . For H (cid:29) a , where a ∼ (cid:112) κ /Y is the effec-tive thickness of the membrane, the elastic energy dom-inates over bending energy, E el (cid:28) E bend . We note that typically, the effective thickness is smaller than the lat-tice spacing, e.g. for graphene a ∼ A . Thus, by ne-glecting E bend and minimizing E el + E b + E vdW (with E b ∼ − P HR ) over H and R , we find H = c R (cid:0) γ/Y (cid:1) / , P = c (cid:0) γ Y (cid:1) / /R. (6)Here the coefficients c ≈ .
86 and c ≈ .
84 has beenobtained from the approximate solution of the Euler-Lagrange equations for h ( r ) and u ( r ) [12]. The re-sults (6) are applicable under conditions γ (cid:28) Y and R (cid:29) a ( Y /γ ) / which guarantee H (cid:28) R and H (cid:29) a ,respectively. We note that the minimization of the en-ergy implies that | E b | ∼ E vdW which gives P ∼ γ/H .This relation between P and H will hold for all consid-ered regimes below. Therefore, in what follows we shallpresent expressions for the maximal height H only. The effect of thermal fluctuations. — A finite tem-perature induces the thermal fluctuations of the mem-brane. These thermal fluctuations are essentially the in-plane and flexural (out-of-plane) phonons. The in-planephonons induce the long-ranged interaction between flex-ural phonons, Eq. (5). The most dangerous are the out-of-plane phonons with wave vectors q < /R ∗ [14, 15],where R ∗ ∼ κ / √ Y T is the so-called Ginzburg length[18]. Therefore, at finite temperature for the bubble ofradius R > R ∗ one needs to integrate out the flexuralphonons with momenta 1 /R < q < /R ∗ before deriva-tion of the Euler-Lagrange equation for h . Essentially,integration over the out-of-plane phonons leads to thesame form of the bending and elastic energies as givenby Eqs. (4) and (5) but with the renormalized bendingrigidity and Young’s modulus [15]: κ ( R ) = κ (cid:0) R/R ∗ (cid:1) η , Y ( R ) = Y (cid:0) R/R ∗ (cid:1) − η . (7)Here η is the universal exponent which depends onthe dimensionality of a membrane and of an embeddedspace. For the clean 2D crystalline membrane in three-dimensional space numerics predicts η ≈ . σ affects the ther-mal fluctuations. There is the characteristic tension σ ∗ = κ /R ∗ ∼ T Y / κ [24, 32–34]. For σ (cid:28) σ ∗ thescaling (7) holds for the interval R ∗ (cid:28) R (cid:28) R σ where R σ = R ∗ ( σ/σ ∗ ) / (2 − η ) is the solution of the equation σ = κ ( R σ ) /R σ . For R (cid:29) R σ the bending rigidity andthe Young’s modulus saturates at the values κ ( R σ ) and Y ( R σ ), respectively.[35] For σ > σ ∗ ( R σ < R ∗ ) the ther-mal fluctuations are completely suppressed and at finitetemperature one can minimize the unrenormalized bend-ing, Eq. (4), and elastic, Eq. (5), energies.The pressure P inside the bubble results in a non-zerotension σ P ∼ P R where R ∼ R /H is the radius of thecurvature of the membrane on the bubble [36]. UsingEq. (6), we find σ P ∼ √ γY . Such tension is enough tosuppress the thermal fluctuations provided σ P (cid:29) σ ∗ , i.e.the standard approach is only valid at low enough tem-peratures: T (cid:28) T γ ∼ κ (cid:112) γ/Y . The energy scale T γ has clear physical meaning of the temperature at whichthe van der Waals energy for the bubble of radius R ∗ be-comes of the order of the temperature. In other words,for T (cid:29) T γ the van der Waals energy does not suppressthe thermal fluctuations. Below we shall study the hightemperature regime, T (cid:29) T γ . Bending dominated regime. — We start from the bubblewith the radius R (cid:28) R ∗ . At such small lengthscale thereis no renormalization of the bending rigidity and Young’smodulus. However, as follows from above, the standardapproach cannot be correct. The only resolution is theassumption that the bending energy is dominated overelastic one, E bend (cid:29) E el , i.e. H (cid:28) a . After minimizationof E bend + E b + E vdW over H and R , we find H = c a (cid:0) T γ /T (cid:1)(cid:0) R/R ∗ (cid:1) , R (cid:28) R ∗ . (8)Here c ≈ .
65 is found from solution of the Euler-Lagrange equation for h ( r ) (Supplemental Material [37]).Now we assume that the bubble radius R (cid:29) R ∗ . Atsuch lengthscales one has to take into account the renor-malization of the bending rigidity and Young’s modulus(provided the scale R σ is large enough). The renormal-ization changes the estimates for the bending and elasticenergies: E bend ∼ κ ( R ) H /R and E el ∼ Y ( R ) H /R .Again we assume that the bending energy is larger thanthe elastic one, E bend (cid:29) E el . This implies that H (cid:28) a ( R/R ∗ ) − η/ . Minimization of E bend + E b + E vdW yields H = c a (cid:0) T γ /T (cid:1)(cid:0) R/R ∗ (cid:1) − η , R ∗ (cid:28) R (cid:28) R ∗ T /T γ , (9)where c ≈ .
90 [37]. The upper bound on R in Eq. (9)comes from the condition E bend (cid:29) E el . In the aboveanalysis we did not take into account the tension of themembrane due to the pressure. Using Eq. (9), we findthe following estimate: σ P ∼ γ ( R/H ) ∼ σ ∗ ( R ∗ /R ) − η ,i.e. R σ ∼ R . Since the power-law renormalization (7)is caused by the flexural phonons with momentum q > /R the tension σ P is indeed irrelevant for the thermalfluctuations in the regime R ∗ (cid:28) R (cid:28) R ∗ T /T γ . Tension dominated regime. — For the bubbles of ra-dius R (cid:29) R ∗ T /T γ the elastic energy is dominated overthe bending one, E el (cid:29) E bend . Then the minimiza-tion of E el + E b + E vdW over H and R implies that E el ∼ | E b | ∼ E vdW . Therefore, the pressure-inducedtension σ P ∼ | E b | /H ∼ E el /H (cid:29) E bend /H . This es-timate means that the pressure induced tension is impor-tant and the corresponding length scale is short, R σ (cid:28) R .In such regime the bending rigidity and Young’s modu-lus are independent of R albeit strongly renormalized.Therefore, we can use the results of the standard ap-proach but with the Young’s modulus Y ( R σ ) instead of Y . In particular, the tension induced by the pressureis given as σ P ∼ (cid:112) γY ( R σ ). Hence the length scale R σ satisfies the following equation: κ ( R σ ) /R σ = (cid:112) γY ( R σ ). FIG. 2. The dependence of the aspect ratio on the radius ofthe bubble at high temperatures, T (cid:29) T γ . Its solution yields R σ ∼ R ∗ T /T γ . This justifies that theprofile of the bubbles with R (cid:29) R ∗ T /T γ are governed bypressure induced tension. The characteristic radius R σ has simple physical meaning. The bubble of such radiushas the adhesion energy, πγR σ , equal to T . Using Eq.(6) with the renormalized Young’s modulus, we find H = c a (cid:0) R/R ∗ (cid:1)(cid:0) T /T γ (cid:1) − η/ , R ∗ T /T γ (cid:28) R. (10)We mention that although the aspect ratio H/R of thebubbles with R (cid:29) R ∗ T /T γ is independent of R , it is notthe constant but depends on the temperature. The valueof the aspect ratio is much larger than one would predicton the basis of the standard approach. The behavior ofthe aspect ratio on R at T (cid:29) T γ is shown in Fig. 2. The anomalous thermodynamics. — The temperaturedependence of the maximal height H depends on theequation of state of the substance inside the bubble. Westart from the case of the liquid bubble. Then we canapproximate the equation of state by the constant vol-ume condition: V = const. We assume that the bubblehas large enough volume, V (cid:29) V γ ∼ a κ /T γ . In theopposite case, V (cid:28) V γ , H is smaller than the effectivethickness of the membrane at all temperatures [37].At low temperatures, T (cid:28) T γ , the maximal height isgiven by Eq. (6), i.e. H is independent of T : H ∼ a (cid:0) V /V γ (cid:1) / , T (cid:28) T γ . (11)At T γ (cid:28) T (cid:28) T γ ( V /V γ ) − η the thermal fluctuations areimportant but the physics is dominated by the tensioninduced by the pressure. Using Eq. (10), we find H ∼ a (cid:32) V T − η V γ T − ηγ (cid:33) , T γ (cid:28) T (cid:28) T γ (cid:18) VV γ (cid:19) − η . (12)At high temperatures, T (cid:29) T γ ( V /V γ ) − η , the maximalheight of the bubble is described by the theory of the tension dominated regime bending dominated regime FIG. 3. The dependence of the maximal height of the bubblewith liquid on the temperature for the case of large volume, V (cid:29) V γ . bending dominated regime, Eq. (9). Then, we find that H is decreasing with increase of temperature: H ∼ a (cid:32) V − η T ηγ V − ηγ T η (cid:33) − η , T γ (cid:18) VV γ (cid:19) − η (cid:28) T. (13)Therefore, in the regime of large volumes, V (cid:29) V γ , themaximal height of the bubble has non-monotonous de-pendence on temperature with the maximum at temper-ature T max ∼ T γ ( V /V γ ) / (4 − η ) (see Fig. 3). The non-monotonous dependence of H implies the change of thesign of the linear thermal expansion coefficient α H attemperature T max : α H = 1 T , T (cid:28) T γ , − η , T γ (cid:28) T (cid:28) T max , − η − η , T max (cid:28) T. (14)Therefore, by measuring the slope of α H against 1 /T onecan extract the bending rigidity exponent of the mem-brane. The result (14) is derived with the neglect oftemperature dependence of the adhesion energy.Now we discuss the case of the bubble with a gas inside.For a sake of simplicity, we use the equation of state of theideal gas, P V = N T , where N is the number of atoms ofthe gas. Using the relations V ∼ HR ∼ γR /P , we findthat the radius of the bubble with the ideal gas is alwaysgiven as R ∼ (cid:112) N T /γ . At low temperature T (cid:28) T γ ,using Eq. (6), we find that the maximal height of thebubble grows with temperature as H ∼ a √ N (cid:0) T /T γ (cid:1) / , T (cid:28) T γ . (15)We note that, strictly speaking, this estimate is valid for T (cid:29) T γ /N . Under assumption of the macroscopic num-ber of atoms, N (cid:29)
1, inside the bubble this limitationon T is completely irrelevant.For high temperatures, T (cid:29) T γ , the bubbles of radius R (cid:29) R ∗ T /T γ (this condition is equivalent to the condi-tion N (cid:29)
1) can be formed only. Therefore, the bubble tension dominated regime
FIG. 4. The dependence of the maximal height of the bubblewith the ideal gas inside on the temperature. is in the tension dominated regime such that its maximalheight is described by Eq. (10). Then, we obtain H ∼ a √ N (cid:0) T /T γ (cid:1) − η/ , T (cid:29) T γ . (16)The above results show that the maximal height of thebubble with the ideal gas inside is the monotonouslygrowing function of temperature (see Fig. 4). There-fore, the linear thermal expansion coefficient is alwayspositive: α H = 1 T (cid:40) / , T (cid:28) T γ , − η/ , T γ (cid:28) T. (17)As in the case of the bubble with a liquid the slope of α H with 1 /T allows to extract the value of the exponent η . Discussion and summary. — With the use of a more real-istic equation of state one can compute the temperaturedependence of the maximal height along the liquid-to-gasisotherm. In particular, one can analyze how the anoma-lous elasticity affects the liquid-to-gas transition in thebubble. This phenomenon has been studied recently inRef. [38] but within the standard approach which ignoresthe thermal fluctuations of the membrane.It is known [36, 39] that the anomalous elasticity ofmembrane affects the stability of spherical membraneshells. The renormalization of the bending rigidity andYoung’s modulus decreases the pressure induced tensiontowards a negative value which is enough for developingthe buckling instability. In principle, similar mechanismof instability is also applicable for the curved membraneabove the bubble. However, our estimates indicate thatsuch buckling instability is out of reach [37].Also it is worthwhile to mention that in the bendingdominated regime there are large thermodynamics fluc-tuations of the bubble height which might complicate theexperimental observation of the predicted dependence ofthe average height of the bubble on R and T [37].Let us estimate the relevant parameters for our theoryin the case of a van der Waals heterostructure made of agraphene monolayer on a monolayer of hBN. Using theknown values of Young’s modulus, bending rigidity, andthe effective thickness of the graphene: Y ≈
22 eV ˚A − , κ ≈ . a ≈ . R ∗ ≈ T = 300 K [32]. Assuming thatthe total adhesive energy is dominated by the adhesiveenergy between graphene and hBN, γ ≈ γ ms ≈ .
008 eV˚A − [2], we find T γ ≈
250 K. This estimate indicatesthat the aspect ratio of the bubbles between grapheneand hBN measured recently [4] can be described by ourtheory in the tension dominated regime, Eq. (10). Usingthe experimentally observed aspect ratio 0 .
11 we obtain T γ ≈
220 K and the adhesion energy γ ≈ .
007 eV ˚A − .The later is 20 per cent higher than the value extracted inRef. [4] within the standard approach, Eq. (6). This im-plies that the proper account for the thermal fluctuationscan be crucial for the precision measurements of the ad-hesion energy via the pressurised blister test. The char-acteristic volume for the bubble between graphene andhBN can be estimated as V γ ≈
10 ˚A . Such smallness ofthe value of V γ suggests that the non-monotonous behav-ior of the maximal height on temperature in graphene-on-hBN structure could be observed experimentally forthe liquid bubbles of few nanometer radius only.In summary, we have demonstrated that the anoma-lous elasticity of membranes affects the profile of a bubblein van der Waals heterostructures at high temperatures.We have shown that the renormalization of the bend-ing rigidity and Young’s modulus results in the anoma-lous dependence of the maximal height of the bubble onits radius and temperature. Our estimates suggest thatfor graphene-based van der Waals heterostructures theanomalous regime predicted in the paper is experimen-tally accessible at ambient conditions. Acknowledgements. —A.A.L. and S.S.B. are gratefulto the “FIZIKA” foundation for support during the 2019Landau Institute summer school on theoretical physicswhere this project was initiated. I.S.B. is grateful to I.Gornyi, V. Kachorovskii, and A. Mirlin for very usefuldiscussions as well as to I. Kolokolov and K. Tikhonovfor valuable remarks. The work was partially supportedby the programs of the Russian Ministry of Science andHigher Education and by the Basic Research Program ofHSE. [1] Z. Dai, L. Liu, and Z. Zhang, “Strain engineering of2d materials: Issues and opportunities at the interface,”Adv. Mater. , 1805417 (2019).[2] Y. T. Megra and J. W. Suk, “Adhesion properties of 2dmaterials,” J. Phys. D: Appl. Phys. , 364002 (2019).[3] S. J. Haigh, A. Gholinia, R. Jalil, S. Romani, L. Britnell,D. C. Elias, K. S. Novoselov, L. A. Ponomarenko, A. K.Geim, and R. Gorbachev, “Cross-sectional imaging ofindividual layers and buried interfaces of graphene-basedheterostructures and superlattices,” Nat. Mater. , 764 (2012).[4] E. Khestanova, F. Guinea, L. Fumagalli, A. K. Geim,and I. V. Grigorieva, “Universal shape and pressure in-side bubbles appearing in van der waals heterostruc-tures,” Nat. Commun. , 12587 (2016).[5] Z. Dai, Y. Hou, D. A. Sanchez, G. Wang, C. J. Brennan,Z. Zhang, L. Liu, and N. Lu, “Interface-governed defor-mation of nanobubbles and nanotents formed by two-dimensional materials,” Phys. Rev. Lett. , 266101(2018).[6] S.M. Ghodsi, C. M. Megaridis, R. Shahbazian-Yassar,and T. Shokuhfar, “Advances in graphene-based liquidcell electron microscopy: Working principles, opportuni-ties, and challenges,” Small Methods , 1900026 (2019).[7] A. V. Tyurnina, D. A. Bandurin, V. G. Khestanova,E.and Kravets, M. Koperski, F. Guinea, A. N. Grig-orenko, A. K. Geim, and I. V. Grigorieva, “Strainedbubbles in van der waals heterostructures as local emit-ters of photoluminescence with adjustable wavelength,”ACS Photonics , 516 (2019).[8] S. P. Koenig, N. G. Boddeti, M. L. Dunn, and J. S.Bunch, “Ultrastrong adhesion of graphene membranes,”Nat. Nanotechnol. , 543 (2011).[9] Narasimha G. Boddeti, Xinghui Liu, Rong Long, Jian-liang Xiao, J. Scott Bunch, and Martin L. Dunn,“Graphene blisters with switchable shapes controlled bypressure and adhesion,” Nano Letters , 6216–6221(2013).[10] D. Lloyd, X. Liu, N. Boddeti, L. Cantley, R. Long, M. L.Dunn, and J. S. Bunch, “Adhesion, stiffness, and insta-bility in atomically thin MoS bubbles,” Nano Letters , 5329 (2017).[11] G. Wang, Z. Dai, J. Xiao, S. Feng, C. Weng, L. Liu,Z. Xu, R. Huang, and Z. Zhang, “Bending of multilayervan der waals materials,” Phys. Rev. Lett. , 116101(2019).[12] K. Yue, W. Gao, R. Huang, and K. M. Liechti, “Ana-lytical methods for the mechanics of graphene bubbles,”J. Appl. Phys. , 083512 (2012).[13] P. Wang, W. Gao, Z. Cao, K. M. Liechti, and R. Huang,“Numerical analysis of circular graphene bubbles,” J.Appl. Mechanics , 040905 (2013).[14] D.R. Nelson and L. Peliti, “Fluctuations in membraneswith crystalline and hexatic order,” Journal de Physique , 1085–1092 (1987).[15] J. A. Aronovitz, “Fluctuations of solid membranes,”Phys. Rev. Lett. , 2634 (1988).[16] M. Paczuski, M. Kardar, and D. R. Nelson, “Landautheory of the crumpling transition,” Phys. Rev. Lett. ,2638 (1988).[17] F. David and E. Guitter, “Crumpling transition in elas-tic membranes: Renormalization group treatment,” Eu-rophysics Lett. (EPL) , 709 (1988).[18] J. Aronovitz, L. Golubovic, and T. C. Lubensky, “Fluc-tuations and lower critical dimensions of crystalline mem-branes,” J. de Physique , 609 (1989).[19] E. Guitter, F. David, S. Leibler, and L. Peliti, “Thermo-dynamical behavior of polymerized membranes,” J. dePhysique , 1787 (1989).[20] P. Le Doussal and L. Radzihovsky, “Self-consistent the-ory of polymerized membranes,” Phys. Rev. Lett. ,1209 (1992).[21] D. Nelson, T. Piran, and S. Weinberg, eds., StatisticalMechanics of Membranes and Surfaces (World Scientific,
Singapore, 1989).[22] P. Le Doussal and L. Radzihovsky, “Anomalous elasticity,fluctuations and disorder in elastic membranes,” Ann.Phys. (N.Y.) , 340 (2018).[23] R. J. T. Nicholl, H. J. Conley, N. V. Lavrik, I. Vlassiouk,Y. S. Puzyrev, V. P. Sreenivas, S. T. Pantelides, andK. I. Bolotin, “The effect of intrinsic crumpling on themechanics of free-standing graphene,” Nat. Commun. ,9789 (2015).[24] I. V. Gornyi, V. Yu. Kachorovskii, and A. D. Mirlin,“Anomalous hooke’s law in disordered graphene,” 2DMaterials , 011003 (2016).[25] H. Ghorbanfekr-Kalashami, K. S. Vasu, R. R. Nair,Fran¸cois M. Peeters, and M. Neek-Amal, “Dependenceof the shape of graphene nanobubbles on trapped sub-stance,” Nat. Commun. (2017).[26] K. Zhang and M. Arroyo, “Coexistence of wrinklesand blisters in supported graphene,” Extreme Mechan-ics Lett. , 23 (2017).[27] M.R. Delfani, “Nonlinear elasticity of monolayer hexag-onal crystals: Theory and application to circular bulgetest,” European Journal of Mechanics - A/Solids , 117(2018).[28] D. A. Sanchez, Z. Dai, A. Wang, Z.and Cantu-Chavez,C. J. Brennan, R. Huang, and N. Lu, “Mechanicsof spontaneously formed nanoblisters trapped by trans-ferred 2D crystals,” PNAS , 7884 (2018).[29] L.D. Landau and E.M. Lifshitz, Theory of Elasticity (Butterworth-Heinemann, 2012).[30] M. J. Bowick, S. M. Catterall, M. Falcioni, G. Thorleif-sson, and K. N. Anagnostopoulos, “The flat phase ofcrystalline membranes,” J. de Physique I , 1321 (1996).[31] A. Tr¨oster, “Fourier monte carlo simulation of crystallinemembranes in the flat phase,” J. Phys.: Conf. Series ,012032 (2013).[32] R. Rold´an, A. Fasolino, K. V. Zakharchenko, andM. I. Katsnelson, “Suppression of anharmonicities incrystalline membranes by external strain,” Phys. Rev.B , 174104 (2011).[33] A. Koˇsmrlj and D. R. Nelson, “Response of thermalizedribbons to pulling and bending,” Phys. Rev. B , 125431(2016).[34] I. S. Burmistrov, I. V. Gornyi, V. Yu. Kachorovskii,M. I. Katsnelson, and A. D. Mirlin, “Quantum elasticityof graphene: Thermal expansion coefficient and specificheat,” Phys. Rev. B , 195430 (2016).[35] Here we neglect weak logarithmic dependence on R of thebending rigidity for R (cid:29) R σ (cf. Ref. [33]). Such logarith-mic corrections are beyond accuracy of our estimates.[36] J. Paulose, G. A. Vliegenthart, G. Gompper, and D. R.Nelson, “Fluctuating shells under pressure,” PNAS ,19551 (2012).[37] See Supplemental Material.[38] P. Zhilyaev, E. Iakovlev, and I. Akhatov, “Liquid–gasphase transition of ar inside graphene nanobubbles on thegraphite substrate,” Nanotechnology , 215701 (2019).[39] A. Koˇsmrlj and D. R. Nelson, “Statistical mechanics ofthin spherical shells,” Phys. Rev. X , 011002 (2017). upplemental Material for “The effect of anomalous elasticity on the profile ofbubbles in van der Waals heterostructures” A. A. Lyublinskaya, S. S. Babkin, and I. S. Burmistrov
2, 3 Moscow Institute of Physics and Technology, 141700, Dolgoprudnyi, Moscow Region, Russia L. D. Landau Institute for Theoretical Physics, Semenova 1-a, 142432, Chernogolovka, Russia Laboratory for Condensed Matter Physics, National Research University Higher School of Economics, 101000 Moscow, Russia
This material contains (i) the accurate analytical calculation of the maximal height of the bubblein the bending dominating regime; (ii) the small volume regime of the liquid bubble; and (iii) theestimates demonstrating the absence of buckling transition of the bubble; (iii) the estimate forthermodynamic fluctuations of the height of the bubble.
S.I. THE MAXIMAL HEIGHT OF THE BUBBLE IN THE BENDING DOMINATING REGIME
Let us introduce the normalized eigenfunctions of the Laplace operator on the disk r R satisfying zero boundaryconditions at r = R : ∆ φ n ( r ) = − ζ n R φ n ( r ) , φ n ( r ) = 1 √ πR | J ( ζ n ) | J (cid:18) ζ n rR (cid:19) . (S1)Here ζ n with n = 1 , , . . . are the zeroes of the zeroth order Bessel function J ( x ). Then, one can write the renormalizedbending energy as follows E bend = 12 Z d r d r ∆ h ( r ) ˆ κ ( r − r )∆ h ( r ) , (S2)where the integral operator ˆ κ is defined via its action on the eigenfunctions of the Laplace operator: Z d r ˆ κ ( r − r ) φ n ( r ) = ζ − ηn κ ( R ) φ n ( r ) . (S3)Here we remind κ ( R ) = κ ( R/R ∗ ) η .Now let us expand the height of the bubble into the series: h ( r ) = P n α n φ n ( r ). This expansion automaticallysatisfies the boundary condition h ( r = R ) = 0. The knowledge of the coefficients α n allows one to find the maximalheight of the bubble: H = P n α n / [ √ πR | J ( ζ n ) | ]. In terms of the coefficients α n the total energy acquires the followingform: E bend + E b + E vdW = κ ( R )2 R X n ζ − ηn α n − √ πP R X n sgn [ J ( ζ n )] ζ n α n + πγR . (S4)Its minimization over α n yeilds α n = 2 √ π P R κ ( R ) sgn [ J ( ζ n )] ζ − ηn . (S5)Hence, for the total energy we obtain E bend + E b + E vdW = − π P R κ ( R ) X n ζ − ηn + πγR . (S6)Minimization with respect to R determines the pressure inside the bubble: P = c p γ κ ( R ) R , c = " (6 − η ) X n ζ − ηn − / ≈ .
72 (S7)Using this result for the pressure, we find the final expression for the maximal height: H = 2 P R κ ( R ) X n ζ − ηn J ( ζ n ) = c R (cid:18) γ κ ( R ) (cid:19) / , c = 2 c X n ζ − ηn J ( ζ n ) ≈ . . (S8)The constants c and c relevant for the low temperature regime R (cid:28) R ∗ are obtained from the results above bysetting η = 0. Then one finds c ≈ .
65 and c ≈ . S.II. THE LIQUID BUBBLE OF A SMALL VOLUME V (cid:28) V γ In the case of the bubble with small volume of liquid, V (cid:28) V γ the temperature dependence of the maximal heightis determined by the bending dominated regime. At low temperatures from Eq. 6 of the main text we obtain thatthe maximal height of the bubble is independent of temperature: H ∼ a (cid:0) V /V γ (cid:1) / , T (cid:28) T γ (cid:0) V /V γ (cid:1) − / . (S9)At high temperatures, as it follows from Eq. 10 of the main text, the maximal height of the bubble starts to decreasewith the temperature: H ∼ a (cid:0) V /V γ (cid:1) − η − η (cid:0) T /T γ (cid:1) − η − η , T γ (cid:0) V /V γ (cid:1) − (cid:28) T. (S10)We note that the power law decay is controlled by the bending rigidity exponent η . The decay of H with increase of T implies the negative linear thermal expansion coefficient α H = (1 /H ) dH/dT . Although interesting on its own, thebubbles with a liquid of small volume V (cid:28) V γ can hardly be detected since as it follows from Eqs. S9 and S10 themaximal height of the bubble is smaller than a . S.III. ABSENCE OF THE BUCKLING INSTABILITY OF A SPHERICAL BUBBLE
In the presence of non-zero curvature radius R of the membrane flexural phonons acquire a true mass equal to Y /R [S1]. This modification of the spectrum of flexural phonons results in negative renormalization of the tension.Perturbatively, a change of the tension δσ can be estimated as [S1]: δσσ ∼ − T Y σR Z d k (2 π ) k ( κ k + σk + Y /R ) . (S11)We emphasize that here the integral over momentum k has infrared cut off due to the radius of the bubble R : k & /R .Eq. S11 can be directly applied to the regime of low temperatures T (cid:28) T γ . Then, using Eq. 6 of the main text,we find δσσ ∼ − TγR ln R √ γY κ . (S12)Since the low temperature theory is applicable for R (cid:29) a ( Y /γ ) / , we find that | δσ | /σ (cid:28) ( T /T γ ) (cid:28) T (cid:29) T γ . For R (cid:28) R ∗ we can use Eq. S11. Then, using Eq. 9of the main text, we obtain ( x = kR ) δσσ ∼ − (cid:18) T γ T (cid:19) (cid:18) RR ∗ (cid:19) Z ∼ dx x ( x + x + γY R / κ ) . (S13)Since the parameter γY R / κ ) ∼ ( T γ /T ) ( R/R ∗ ) (cid:28)
1, integrating over x , we find δσσ ∼ − (cid:18) T γ T (cid:19) (cid:18) RR ∗ (cid:19) . (S14)Therefore we find that | δσ | /σ (cid:28) ( T γ /T ) (cid:28) R (cid:28) R ∗ . In the case of intermediate radius R ∗ (cid:28) R (cid:28) R σ = R ∗ T /T γ , one needs to modify Eq. S11 in order to include renormalization of the bending rigidity and Young’s modulus: δσσ ∼ − T Y ( R ) σR Z d k (2 π ) k ( κ ( R ) k + σk + Y ( R ) /R ) . (S15)Using Eq. 10 of the main text, we obtain δσσ ∼ − (cid:18) RR σ (cid:19) Z ∼ dx x ( x + x + γY ( R ) R / κ ( R )) . (S16)Since γY ( R ) R / κ ( R ) ∼ ( R/R σ ) (cid:28)
1, we obtain | δσ | σ ∼ (cid:18) RR σ (cid:19) (cid:28) , R ∗ (cid:28) R (cid:28) R σ . (S17)Finally, we consider the case R (cid:29) R σ . Here we can use Eq. S12 with κ and Y substituted by κ ( R σ ) and Y ( R σ ),respectively. Then, we find | δσ | σ ∼ (cid:18) R σ R (cid:19) ln RR σ (cid:28) , R (cid:29) R σ . (S18)The above results demonstrates that for both low and high temperature regimes, the renormalization of the tensiondue to the finite curvature is weak and cannot lead to the buckling instability. The only exception is the bubbles withthe radius R ∼ R σ for which our estimates give | δσ | /σ ∼
1. This implies that for R ∼ R σ the buckling instabilitycould occur in principle. In order to resolve this issue one needs to perform more accurate RG scheme similar tothe one of Ref. [S2]. However, since for R (cid:28) R σ and R (cid:29) R σ the buckling instability is absent the scenario withinstability for R ∼ R σ seems to be unlikely. S.IV. THE THERMODYNAMIC FLUCTUATIONS OF THE HEIGHT OF THE BUBBLE
Let us estimate the thermodynamic fluctuations of the height of the bubble in the bending dominated regime, R (cid:28) R ∗ T /T γ . Using Eq. 4 of the main text, we find h ( δH ) i ∼ T Z d k (2 π ) κ ( R ) k . (S19)As above the integral over momentum k has infrared cut off due to the radius of the bubble R : k & /R . Hence wefind that the dispersion of height fluctuations is given as follows h ( δH ) i ∼ T R / κ ( R ) . (S20)Then for R (cid:28) R ∗ T /T γ we obtain the following estimate: p h ( δH ) i H ∼ TT γ R ∗ R (cid:29) . (S21)This inequality implies that there are large thermodynamic fluctuations of the height of the bubble in the bendingdominated regime.For the tension dominated regime R (cid:29) R σ = R ∗ T /T γ one needs to take into account the tension induced by thepressure: h ( δH ) i ∼ T Z d k (2 π ) κ ( R σ ) k + σ P k ∼ T R σ κ ( R σ ) ln RR σ . (S22)Then we obtain the following estimate: p h ( δH ) i H ∼ TT γ R ∗ R ln RR σ (cid:28) . (S23)Therefore the thermodynamic fluctuations of the height of the bubble in the tension dominated regime are stronglysuppressed. [S1] J. Paulose, G. A. Vliegenthart, G. Gompper, and D. R. Nelson, Fluctuating shells under pressure , PNAS , 19551(2012).[S2] A. Kosmrlj and D. R. Nelson,
Statistical mechanics of thin spherical shells , Phys. Rev. X7