The effect of anomalous elasticity on the bubbles in van der Waals heterostructures
A. A. Lyublinskaya S. S. Babkin I. S. Burmistrov






 

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 a bubble in van der Waals heterostructures. Our theory generalizes the non-linear plate theory as well as membrane theory of the pressurised blister test to incorporate the power-law scale dependence of the bending rigidity and Young's modulus of a two-dimensional crystalline membrane. This scale dependence 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 this anomalous elasticity affects dependence of the maximal height of the bubble on its radius and temperature. We identify the characteristic temperature above which the anomalous elasticity is important. It is suggested that for graphene-based van der Waals heterostructures the predicted anomalous regime is experimentally accessible at the room temperature.


The effect of anomalous elasticity on the bubbles in van der Waals heterostructures
    I. S. Burmistrov
    December 30, 2023
===================================================================================



 — Mechanical properties of two-dimensional  (2D) materials, especially, of van der Waals heterostructures, have recently attracted a lot of interest  in view of their potential applications <cit.>. The simplest example of van der Waals heterostructure is two monolayers, e.g.
 graphene, hexagonal boron nitride (hBN), MoS_2, assembled together. Strong adhesion between monolayers <cit.> results in atomically clean interfaces in which all contaminating substances are combined into bubbles <cit.>. Recently, these bubbles inside van der Waals heterostructures have been studied experimentally <cit.>. Similar bubbles are formed between an atomic monolayer and a substrate, e.g. SiO_2 <cit.>. There are many suggestions of  practical usage of the bubbles inside van der Waals heterostructures, for example, the 
graphene liquid cell microscopy <cit.>, 
controlled room-temperature photoluminescence emitters <cit.>, etc.


 
The mechanics of monolayers due to these bubbles is considered to be analogous to the one of the pressurized blister test  which has recently become the routine method to measure simultaneously the Young's modulus and adhesion energy of a monolayer on a substrate <cit.>. Usually, the pressurized blister test is described either by  nonlinear plate model or by membrane theory (see e.g. <cit.>). These standard elastic theories of deformed plates ignore the fact that elastic properties of an atomic monolayer are those of 2D crystalline membranes <cit.>. The most striking feature of mechanics of membranes is anomalous elasticity which results in non-linear (so-called, anomalous) Hooke's law for small applied stress (see Refs. <cit.> for a review).
Recently, this anomalous Hooke's law has been measured in graphene <cit.>. However, until present, see e.g. Refs. <cit.>, the anomalous elasticity of membranes is completely ignored in description of mechanical properties of monolayers in the presence of the bubbles.

The present paper generalizes the classical theory of the pressurised blister test to incorporate the anomalous elasticity of a membrane. Our approach explicitly takes into account the power-law renormalization of the bending rigidity and Young's modulus. It is shown that above a certain temperature the dependence 
of the bending rigidity and Young's modulus of a membrane on the radius  of the bubble results in the non-analytic dependence of its maximal height on the radius and temperature.  




 — The model for the description of the profile of a bubble between the membrane  and the substrate (see Fig. <ref>) is well established <cit.>. It can be formulated in terms of the (free) energy which is the sum of the following four terms:    

    E=E_ bend+E_ el+E_ b+E_ vdW .

Here the first contribution describes the energy cost related with bending of the membrane:

    E_ bend = ϰ_0/2∫ d^2 r[ (Δ h)^2 + (Δu)^2 ] ,

where ϰ_0 denotes the bare bending rigidity of the membrane.  Here u={u_x,u_y} and h are the in-plane and out-of-plane displacements of the membrane (see Fig. <ref>). The second term is the standard elastic energy <cit.>:

    E_ el = ∫ d^2 r (
    μ_0 u_αβu_βα
    +λ_0 u_ααu_ββ /2  ) ,

where μ_0 and λ_0 stand for the Lamé coefficients and  u_αβ =  ( ∂_β u_α + ∂_α u_β  + ∂_α h ∂_β h + ∂_αu∂_βu )/2 is the strain tensor.
The  third term E_ b describes the contribution of the bubble substance. Under assumption of the constant pressure P inside the bubble one has E_ b = - P V, where the bubble volume can be approximated as 
V = ∫ d^2 r   h(r). The last term in Eq. (<ref>) describes the van der Waals interaction between the membrane and the substrate in the presence of the bubble. It can be approximated as
E_ vdW = πγ R^2.  Here R denotes the radius of the bubble (see Fig. <ref>) and  γ= γ_ ms-γ_ mb-γ_ sb where γ_ ms, γ_ mb, and γ_ sb
are the adhesion energies between the membrane and substrate, between the membrane and the substance inside the bubble, and between the substrate and the substance, respectively. The form (<ref>) of the total energy is well justified if the maximal height of the bubble, H=h(0), is small compared to the radius, 

H ≪ R.


Throughout the paper we shall assume that this condition is fulfilled.



 — In order to compute the profile of the spherical bubble, initially, one needs to solve the Euler-Lagrange equations for u(r) and h(r) with the proper boundary conditions and compute the energy E as a function of R and H. Usually, instead of the solution of the Euler-Lagrange equations the approximate solutions either within the non-linear plate theory or within the membrane theory are used, see e.g. Ref. <cit.>. Finally, one have to minimize E with respect to the both R and H. The minimization procedure allows one to find the maximal height H and the pressure P as a function of 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 horizontal deformation: u_ max∼ H^2/R. In the considered regime, H/R ≪ 1,

 the horizontal displacement is also small, u_ max≪ H. This implies that one can neglect the term ∂_αu∂_βu in u_αβ
. Also this allows one to omit the term (Δu)^2 in the bending energy such that it reads:

    E_ bend = ϰ_0/2∫ d^2 r (Δ h)^2 .


In the absence of ∂_αu∂_βu in u_αβ the elastic energy becomes quadratic in u. This implies that the Euler-Lagrange equation for u(r) become linear and can be solved for arbitrary configuration of h(r) (even not necessarily obeying the Euler-Lagrange equation). In other words, a horizontal deformation is adjusted to any vertical displacement. Therefore, E_ el is given as <cit.>:


    E_ el = Y_0/8∫d^2r [K_αα
    -
    ∂_α∫ d^2r^'𝒢(r,r^')
    ∂_β
    K_αβ(r^')  ]^2 ,

where K_αβ = ∂_α h ∂_β h and Y_0 = 4μ_0(μ_0+λ_0)/2μ_0+λ_0 is the Young's modulus. The function 𝒢(r,r^') is the Green's function of the Laplace operator on the disk r⩽ R.


Using Eqs. (<ref>) and (<ref>), we can estimate the bending and elastic energies as 
E_ bend∼ϰ_0 H^2/R^2 and
E_ el∼ Y_0 H^4/R^2. For H ≫ a, where a ∼√(ϰ_0/Y_0) is the effective thickness of the membrane, the elastic energy dominates over bending energy, E_ el≪ E_ bend. We note that typically, the effective thickness is smaller than the lattice spacing, e.g.  for graphene a ∼ 1 Å. Thus, by neglecting E_ bend and minimizing E_ el+ E_ b + E_ vdW (with E_ b∼ -P H R^2) over H and R, we find

    H = c_1 R (γ/Y_0 )^1/4,     
    P = c_2 (γ^3 Y_0)^1/4/R .

Here the coefficients c_1 ≈ 0.86 and c_2≈ 1.84 has been obtained from the approximate solution of the Euler-Lagrange equations for h(r) and u(r) <cit.>.
The results (<ref>) are applicable under conditions γ≪ Y_0 and R≫ a (Y_0/γ)^1/4 which guarantee H ≪ R and H≫ a, respectively.
We note that the minimization of the energy implies that |E_ b| ∼ E_ vdW which gives P∼γ/H. This relation between P and H will hold for all considered regimes below. Therefore, in what follows we shall present expressions for the maximal height H only. 



 — A finite temperature induces the thermal fluctuations of the membrane. These thermal fluctuations are essentially the in-plane and flexural (out-of-plane) phonons. The in-plane phonons  induce the long-ranged interaction between flexural phonons, Eq. (<ref>). The most dangerous are the out-of-plane phonons with wave vectors q< 1/R_* <cit.>, where R_* ∼ϰ_0/√(Y_0 T) is the so-called Ginzburg length <cit.>. Therefore, at finite temperature for the bubble of radius R>R_*  one needs to integrate out the flexural phonons with momenta 1/R<q<1/R_* before derivation of the Euler-Lagrange equation for h. Essentially, integration over the out-of-plane phonons leads to the same form of the bending and elastic energies as given by Eqs. (<ref>) and (<ref>) but with the renormalized bending rigidity and Young's modulus <cit.>:

    ϰ(R) = ϰ_0 (R/R_*)^η,      Y(R)=Y_0 (R/R_*)^-2+2η .

Here η is the universal exponent which depends on the dimensionality of a membrane and of an embedded space. For the clean 2D crystalline membrane in three-dimensional space numerics predicts η≈ 0.8 <cit.>. 

 
The presence of a non-zero tension σ affects the thermal fluctuations. There is the characteristic tension σ_* = ϰ_0/R_*^2 ∼ TY_0/ϰ_0
<cit.>.
For σ≪σ_* the scaling (<ref>) holds for the interval R_*≪ R≪ R_σ where R_σ = R_* (σ/σ_*)^1/(2-η) is the solution of the equation σ = ϰ(R_σ)/R_σ^2. For R≫ R_σ the bending rigidity and the Young's modulus saturates at the values ϰ(R_σ) and Y(R_σ), respectively.[Here we neglect weak logarithmic dependence on R of the bending rigidity for R≫ R_σ (cf. Ref. <cit.>). Such logarithmic corrections are beyond accuracy of our estimates.]  
For σ > σ_* (R_σ<R_*) the thermal fluctuations are completely suppressed and at finite temperature one can minimize the unrenormalized bending, Eq. (<ref>),
and elastic, Eq. (<ref>), energies.

The pressure P inside the bubble results in a non-zero tension σ_P ∼ P R_0 where   
R_0 ∼ R^2/H is the radius of the curvature of the membrane on the bubble <cit.>. Using Eq. (<ref>), we find σ_P ∼√(γ Y_0). Such tension is enough to suppress 
the thermal fluctuations provided σ_P ≫σ_*, i.e. the standard approach is only valid at low enough temperatures:

T ≪  T_γ∼ϰ_0 √(γ/Y_0) .


The energy scale T_γ has clear physical meaning of  the temperature at which the van der Waals energy for the bubble of radius R_* becomes of the order of the temperature. In other words, for T ≫ T_γ the van der Waals energy does not suppress the thermal fluctuations. Below we shall study the high temperature regime, T ≫ T_γ. 



 — We start from the bubble with the radius R≪ R_*. At such small lengthscale there is no renormalization of the bending rigidity and Young's modulus. However, as follows from above, the standard approach cannot be correct. The only resolution is the assumption that the bending energy is dominated over elastic one, E_ bend≫ E_ el, i.e. H ≪ a. After minimization of E_ bend+E_ b+E_ vdW over H and R, we find

    H = c_3 a  (T_γ/T )  (R/R_*)^2 ,
         R≪ R_* .

Here c_3≈ 0.65 is found from solution of the Euler-Lagrange equation for h(r) (Supplemental Material  <cit.>). 

 
Now we assume that the bubble radius R≫ R_*. At such lengthscales one has to take into account the renormalization of the bending rigidity and Young's modulus (provided the scale R_σ is large enough). The renormalization changes the estimates for the bending and elastic energies:
E_ bend∼ϰ(R) H^2/R^2 and E_ el∼ Y(R) H^4/R^2.
Again we assume that the bending energy is larger than the elastic one, E_ bend≫ E_ el. This implies that H≪ a (R/R_*)^1-η/2. Minimization of E_ bend+E_ b+E_ vdW yields 

    H  = c_4 a  (T_γ/T ) (R/R_*)^2-η/2,    R_*≪ R≪ R_* T/T_γ ,

where c_4≈ 0.90 <cit.>. 
The upper bound on R in Eq. (<ref>) comes from the condition  E_ bend≫ E_ el. In the above analysis we did not take into account the tension of the membrane due to the pressure. Using Eq.  (<ref>), we find
the following estimate: σ_P ∼γ (R/H)^2∼σ_* (R_*/R)^2-η, i.e. R_σ∼ R. Since the power-law renormalization (<ref>) is caused by the flexural phonons with momentum q>1/R the tension σ_P is indeed irrelevant for the thermal fluctuations in the regime R_*≪ R ≪ R_* T/T_γ.  



 —  For the bubbles of radius R≫ R_* T/T_γ the elastic energy is dominated over the bending one, 
E_ el≫ E_ bend. Then the minimization of E_ el+E_ b+E_ vdW over H and R implies that E_ el∼ |E_ b| ∼ E_ vdW. Therefore, the pressure-induced tension σ_P ∼ |E_ b|/H^2∼ E_ el/H^2 ≫ E_ bend/H^2. This estimate means that the pressure induced tension is important and the corresponding length scale is short, R_σ≪ R. In such regime the bending rigidity and Young's modulus are independent of R albeit strongly renormalized. Therefore, we can use the results of the standard approach but with the Young's modulus Y(R_σ) instead of Y_0. In particular, the tension induced by the pressure is given as σ_P ∼√(γ Y(R_σ)). Hence the length scale R_σ satisfies the following equation: ϰ(R_σ)/R_σ^2=√(γ Y(R_σ)). Its solution yields
R_σ∼ R_* T/T_γ. This justifies that the profile of the bubbles with R≫ R_* T/T_γ are governed by pressure induced tension. The characteristic radius R_σ has simple physical meaning. The bubble of such radius has the adhesion energy, πγ R_σ^2, equal to T. Using Eq. (<ref>) with the renormalized Young's modulus, we find 

    H = c_1 a  (R/R_* ) (T/T_γ)^-η/2,    R_*T/T_γ≪ R .

We mention that although the aspect ratio H/R of the bubbles with R≫ R_*T/T_γ is independent of R, it is not the constant but depends on the temperature. 
The value of the aspect ratio is much larger than one would predict on the basis of the standard approach. The behavior of the aspect ratio on R at T≫ T_γ is shown in Fig. <ref>. 
 
 
 
 
 
 — The temperature dependence of the maximal height H depends on the equation of state of the substance inside the bubble. 

We start from the case of the liquid  bubble. Then we can approximate the equation of state 
by the constant volume condition: V= const. We assume that the bubble has large enough volume, V≫ V_γ∼ a^3 ϰ_0/T_γ. In the opposite case, V≪ V_γ,  

H is smaller than the effective thickness of the membrane at all temperatures <cit.>.

At low temperatures, T≪ T_γ, the maximal height is given by Eq. (<ref>), i.e. H is independent of T: 

    H ∼ a  (V/V_γ )^1/3,     T
    ≪ T_γ  .
 
At 
T_γ≪ T≪ T_γ (V / V_γ)^2/4-η the thermal fluctuations are important but the physics is dominated by the tension induced by the pressure. Using Eq. (<ref>), we find 


    H ∼ a  (V T^1-η/V_γ T_γ^1-η )^1/3
    ,  
    T_γ≪
    T ≪ T_γ (V/V_γ )^2/4-η .
 
At high temperatures, T≫ T_γ (V / V_γ)^2/4-η, the maximal height of the bubble is described by the theory of the bending dominated regime, Eq. (<ref>). Then, we find that H is decreasing with increase of temperature:

    H ∼ a  (V^4-η T^η_γ/V_γ^4-ηT^η )^1/8-η
    ,     
    T_γ (V/V_γ )^2/4-η≪ T  .
 
Therefore, in the regime of large volumes, V≫ V_γ, the maximal height of the bubble has non-monotonous dependence on temperature with the maximum 
at temperature T_ max∼ T_γ (V/V_γ)^2/(4-η) (see Fig. <ref>). The non-monotonous dependence of H implies the change of the sign of the linear thermal expansion coefficient α_H at temperature T_ max:

    α_H = 1/T
    0,     T≪ T_γ , 
    1-η/3,     T_γ≪ T≪ T_ max, 
    
    -η/8-η,     T_ max≪ T .

Therefore, by measuring the slope of α_H against 1/T one can extract the bending rigidity exponent of the membrane. The result (<ref>) is derived with the neglect of temperature 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 the ideal gas, PV=N T, where N is the number of atoms of the gas. Using the relations V∼ HR^2∼γ R^2/P, we find that the radius of the bubble with the ideal gas is always given as R∼√(N T/γ). At low temperature T≪ T_γ, using  Eq. (<ref>), we find that the maximal height of the bubble grows with temperature as

    H ∼ a √(N) (T/T_γ )^1/2,      T≪ T_γ .
  
We note that, strictly speaking, this estimate is valid for T≫ T_γ/N. Under assumption of the macroscopic number of atoms, N≫ 1, inside the bubble this limitation on T is completely irrelevant.   

For high temperatures, T≫ T_γ, the bubbles of radius R≫ R_* T/T_γ (this condition is equivalent to the condition N≫ 1) can be formed only. Therefore, the bubble is in the tension dominated regime such that its maximal height is described by Eq. (<ref>). Then, we obtain

    H ∼ a √(N) (T/T_γ )^1-η/2,      T≫ T_γ .

The above results show that the maximal height of the bubble with the ideal gas inside is the monotonously growing function of temperature (see Fig. <ref>). Therefore, the linear thermal expansion coefficient is always positive:

    α_H = 1/T
    1/2,         T≪ T_γ , 
    
    1-η/2,         T_γ≪ T .

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

 



 — With the use of a more realistic equation of state one can compute the temperature dependence of the maximal height along the liquid-to-gas isotherm. In particular, one can analyze how the anomalous elasticity affects the liquid-to-gas transition in the bubble. This phenomenon has been studied recently in Ref. <cit.> but within the standard approach which ignores the thermal fluctuations of the membrane.  

It is known <cit.> that the anomalous elasticity of membrane affects the stability of spherical membrane shells. The renormalization of the bending rigidity and Young's modulus decreases the pressure induced tension towards a negative value which is enough for developing the buckling instability. In principle, similar mechanism of instability is also applicable for the curved membrane above the bubble. However, our estimates indicate that such buckling instability is out of reach <cit.>.  

Also it is worthwhile to mention that in the bending dominated regime there are large thermodynamics fluctuations of the bubble height which might complicate the experimental observation of the predicted dependence of the average height of the bubble on R and T <cit.>.

Let us estimate the relevant parameters for our theory in the case of a van der Waals heterostructure made of a graphene monolayer on a monolayer of hBN. Using the known values of Young's modulus, bending rigidity, and the effective thickness of the graphene: Y_0 ≈ 22 eV Å^-2,  
ϰ_0 ≈ 1.1 eV, a≈ 0.6 Å, we can estimate the Ginzburg length as R_* ≈  4 Å  at T=300 K <cit.>. Assuming that 
the total adhesive energy is dominated by the adhesive energy between graphene and hBN, γ≈γ_ ms≈ 0.008 eV Å^-2 <cit.>, we find  T_γ≈ 250 K. This estimate indicates that the aspect ratio of the bubbles between graphene and hBN measured recently <cit.>  can be described by our theory in the tension dominated regime, Eq. (<ref>). Using the experimentally observed aspect ratio 0.11 we obtain T_γ≈ 220 K and the adhesion energy γ≈ 0.007 eV Å^-2. The later is 20 per cent higher than the value extracted in Ref. <cit.> within the standard approach, Eq. (<ref>). This implies that the proper account for the thermal fluctuations can be crucial for the precision measurements of the adhesion energy via the pressurised blister test. The characteristic volume for the bubble between graphene and hBN can be estimated as  V_γ≈ 10 Å^3. Such smallness of the value of V_γ suggests that the non-monotonous behavior of the maximal height on temperature in graphene-on-hBN structure could be observed experimentally for the liquid bubbles of few nanometer radius only. 




In summary, we have demonstrated that the anomalous elasticity of membranes affects the profile of a bubble in van der Waals heterostructures at high temperatures.

We have shown that the renormalization of the bending rigidity and Young's modulus results in the anomalous dependence of the maximal height of the bubble on its radius and temperature. 

Our estimates suggest that for graphene-based van der Waals heterostructures the anomalous regime predicted in the paper is experimentally accessible at ambient conditions. 

 —A.A.L. and S.S.B. are grateful to the “FIZIKA” foundation for support during the 2019 Landau Institute summer school on theoretical physics where this project was initiated. I.S.B. is grateful to I. Gornyi, V. Kachorovskii, and A. Mirlin for very useful discussions as well as to I. Kolokolov and K. Tikhonov for valuable remarks. The work was partially supported by the programs of the Russian Ministry of Science and Higher Education and by the Basic Research Program of HSE.






in 1,...,3 

[pages=,turn=false]supplement.pdf


