The extraordinary mechanical properties of two-dimensional (2D) materials have attracted considerable interest from the physics and material science communities (123–4). Such materials present a very versatile playground allowing realization of multiple modeling objects, such as the thinnest possible membranes with miniscule bending rigidity or very stiff plates, all with very reproducible parameters. Furthermore, both the electronic and mechanical properties of such crystals can be fine-tuned by assembling them into van der Waals (vdW) heterostructures (56–7), which allow realization of yet even broader range of materials parameters and boundary conditions. Of special interest are the properties of free-standing mono- and few-layer 2D crystals: membranes and bubbles (89–10). The elastic response of such systems is highly nontrivial, being governed by the vdW membrane-substrate’s attraction, the bending rigidity of 2D solid membranes, and their Young modulus. Furthermore, atomic thinness of the membranes makes it possible to switch between different regimes—from where the bending rigidity is much smaller than the Young modulus (monolayer membrane), to where they are comparable (in the case of few-layer membranes). A common outcome of this interplay is the instability of the neutrally planar state and the consequent emergence of complex wrinkle patterns when the membrane experiences minute levels of compression (1112–13). Understanding the conditions under which such instabilities emerge and spread throughout the system is of utmost importance both fundamentally and for the functionality of various vdW heterostructure devices. Specifically, the formation of wrinkling patterns is usually associated with a nontrivial strain distribution (14, 15), which, for piezoelectric 2D crystals such as MoS₂ and hBN (1617–18), leads to a complex distribution of the electric field and can be used in respective devices. Furthermore, a careful analysis of the emerging patterns can be harnessed for metrological purposes—extracting valuable information on the material parameters that characterize the elastic moduli of 2D membranes and their vdW interaction with the substrate.

In this work, we study patterns of radially oriented wrinkles induced by bubbles in vdW heterostructures that consist of hexagonal boron nitride (hBN) or molybdenum disulfide (MoS₂) layers on top of another 2D material. We show that the periodicity of the emerging patterns reflects a balance between the membranal bending rigidity (B) and the steepness of the vdW potential well, namely, its second derivative ($$ {V_{vdW}}^{\text{'}\text{'}}$$), evaluated at the equilibrium value of the membrane-substrate distance (19). Although this type of balance between bending rigidity and stiffness of an “effective substrate” has been known to govern wrinkle patterns in thin polymer sheets (14, 20212223–24), the current report provides a noninvasive probe for the vdW interaction between atomically thin membranes and atomically flat substrates. Furthermore, the relative crystallographic orientation between the crystalline, atomically thin membrane and the substrate can be controlled precisely, giving us an extra knob to change the vdW interaction. Our results demonstrate that classical nonlinear mechanics of solids, which describes strain-induced elastic instabilities, can be applied (with some modifications) to the description of pattern formation in complex vdW heterostructures and provide a metrological tool for characterizing the interplay between the vdW attraction and the bending modulus of monolayers and multilayer composites.

We produced high-quality vdW heterostructures with atomically flat monocrystals of hBN on top of graphene (both in incommensurate and commensurate configurations) and MoS₂ on top of hBN using the standard dry transfer technique (25). In brief, we mechanically exfoliate and identify our top 2D material on a polymer film consisting of two sacrificial layers, one of polymethylglutarimide (PMGI) and the other of poly(methyl methacrylate) (PMMA). We develop the PMGI layer from beneath the PMMA layer to create a free-standing membrane that we could easily manipulate with the crystal on top. Then we inverted the membrane and we positioned it above our substrate 2D material, which is resting on a SiO₂/Si substrate, using a set of micromanipulation stages with spatial accuracy better than 5 μm and alignment precision ∼0.5° to control the twist angle. For the commensurate configuration, we aligned the crystal lattices by choosing straight edges of both top and substrate 2D crystals, which indicate the principal crystallographic directions, and brought them into contact making their edges parallel. For the incommensurate configuration, we rotated the top crystal by ∼15° with respect to the 2D substrate edges before bringing them into contact. We then removed the PMMA by simply peeling back the membrane. In such structures, bubbles of trapped substances are spontaneously formed between the substrate and the membrane, such that the membrane portion on top of the bubble is completely delaminated from the substrate. Bubbles appear when 2D crystals with high enough adhesion are brought together due to a self-cleansing mechanism (7, 26), which pushes trapped contaminants to gather into bubbles (9, 27). We imaged our samples with atomic force microscopy (AFM) (28, 29) using a dynamic mode (30, 31), so the images were taken in a gentle noninvasive way ensuring that the sample was not damaged or deformed. We found that the typical radii R of bubbles are in the range 50–300 nm, and their heights H are 10–30 nm; however, the aspect ratio H/R is constant for each given membrane thickness. This constant ratio was shown to be $$ C_0{\left(\frac{\text{Γ}}{Y}\right)}^{1/4}$$, where $$ C_0$$ is a numerical constant, Γ is the membrane-substrate surface energy (i.e., magnitude of vdW attraction), and $$ Y\gg \text{Γ}$$ is the 2D membranal stretching modulus (9).

We start our analysis with the incommensurate configuration. Fig. 1 shows the development of a “corona” of wrinkles (azimuthal undulations of the shape at the laminated portion around the bubble) for membranes of different thicknesses.

From profiles around the topographical images of the wrinkles around the bubbles we can determine the wavelength $$ \ell $$ of the wrinkles as $$ \ell \left(r\right)=2\pi r/\hspace{0.17em}m\left(r\right)$$, where $$ m\left(r\right)$$ is the number of azimuthal undulations at a distance r. Fig. 1 C and D show examples of the wavelength determination for the case of bilayer hBN for two different radii, yielding similar values (see SI Appendix, Fig. S1 for information on how the wrinkles were counted).

The mere existence of radial wrinkles is attributed to the pressure in the bubble. Rather than a uniform, isotropic surface tension (as would be for a gas-filled vesicle), such a pressure gives rise to an anisotropic stress field in the membrane, whose radial component is tensile, $$ {\sigma }_{rr}\sim {\text{Γ}}^{1/2}{Y}^{1/2}\gg \text{Γ}$$, provided that $$ Y\gg \text{Γ}$$. The radial force is pulling inward the portion of the membrane that remains attached to the substrate $$ \left(r>R\right)$$, and consequently causing there azimuthal (hoop) compression, $$ {\sigma }_{\theta \theta }<0$$ (9). Since they emerge just in order to suppress such a compression, the radial extent of wrinkles corresponds to the region of azimuthal compression, which is known to be proportional to the bubble’s radius R, where the proportionality constant is the ratio $$ {\sigma }_{rr}\left(r\right)/{\sigma }_{\infty }$$ between the radial tension $$ {\sigma }_{rr}\left(r\right)\sim {\text{Γ}}^{1/2}{Y}^{1/2}$$ that pulls inward at the bubble’s edge and the far-field (uniform, isotropic) stress $$ {\sigma }_{\infty }$$ in the membrane (24). This is illustrated in Fig. 2, where we show that the formation of these radially oriented wrinkles is correlated with the presence of thin layers of trapped molecules in the vicinity of the bubbles, as observed by monitoring the same bubbles just after fabrication in ambient conditions and after leaving them overnight at high humidity (relative humidity ≳ 90%) (Fig. 2 A–D). A plausible reason for the strong effect of humidity, schematically depicted in Fig. 2 E–H, is that the intake of water molecules between the top and the bottom layers acts to increase the pressure inside the bubble, thereby enhancing the radial tension in the top membrane and, consequently, the hoop compression induced by it, promoting the development of radial wrinkles around the bubble.

Throughout the azimuthally compressed corona of the bubble, the low bending rigidity of the membrane underlies an elastic instability, whereby the membrane seeks to deflect away from the substrate in order to suppress hoop compression and eliminate the large energetic cost it entails. If the vdW attraction were sufficiently weak, one may expect this instability to give rise to a highly nonuniform pattern, where a complete detachment from the substrate in a few narrow angular sectors enables the membrane to retain a full contact with the substrate elsewhere. Such a delamination pattern would be characterized by a small number of narrow, high-amplitude and radially oriented “folds,” separated by perfectly flat, nonundulatory zones, as was assumed in ref. 11, but in sharp contrast with our observations. Instead, the emergence of a simple periodic pattern, characterized by a single wavelength, suggests another instability mechanism, whereby the vdW attraction is sufficiently strong to resist delamination anywhere but beneath the highly pressurized bubble, and small-amplitude radial wrinkles emerge to suppress hoop compression in the corona (see schematic Fig. 3A). This phenomenon resembles the effect of a liquid drop placed on top of a floating polymer sheet, where the Laplace pressure inside the drop gives rise to preferentially radial tension in the sheet, and to the suppression of would-be hoop compression through radial wrinkles (21, 24).

The typical “wavelength” $$ \ell $$ of these wrinkles depends on the thickness of the membranes (see Fig. 4A). Denoting the membrane-substrate distance by $$ z\left(r,\theta \right)$$, the wrinkle pattern can be modeled as:

where $$ A\left(r\right)$$ and $$ m\left(r\right)=2\pi r/\ell \left(r\right)$$ are, respectively, the amplitude and number of azimuthal undulations at a distance r > R from the center of the bubble. From the AFM images we can characterize the shape $$ z\left(r,\theta \right)$$ and observe a few noteworthy features:

1)The number of wrinkles decreases sharply with the number of layers N.

2)The wrinkle number $$ m\left(r\right)$$ increases with radial distance r. Furthermore, a detailed inspection indicates that $$ m\left(r\right)\propto r$$, such that wrinkles are “nucleated” at various locations (32), and the wavelength $$ \ell \left(r\right)=2\pi r/m\left(r\right)$$ is close to a constant value, $$ \ell \left(r\right)={\ell }_0$$ (SI Appendix, Fig. S2).

3)In the vicinity of the bubble, $$ r\hspace{0.17em}\ge R$$, the nonoscillatory component of the shape, $$ z_0\left(r\right)$$, decays exponentially, such that most of the wrinkled zone is nearly flat in the radial direction and highly curved in the azimuthal direction (SI Appendix, Fig. S3).

4)Notwithstanding their high azimuthal curvature, the amplitude of wrinkles, $$ A\left(r\right)$$, is smaller than 1 nm in the majority of the wrinkled zone (SI Appendix, Fig. S3).

5)For a given type (and number of layers) of the top membrane, the average wavelength $$ {\ell }_0$$ does not depend on the size (radius or height) of bubbles (Fig. 4 A–C).

Focusing our attention on the wavelength of wrinkles, we recall a general relation, known as the “local wavelength law” (14, 23):

Here $$ B_{eff}$$ is the effective bending rigidity of the membrane and $$ K_{eff}\left(r\right)$$ is an “effective stiffness.” Eq. 2 states that the wavelength of wrinkles reflects an energetic balance between the bending rigidity of the membrane, which favors undulations with small curvature (hence large $$ \ell $$), and the restoring forces that favor small amplitude (hence small $$ \ell $$). Physically, effective stiffness can originate from multiple contributions, such as attraction to the substrate (vdW interaction), as well as tension- and curvature-induced stiffness (associated with energy cost for deformation of curved membranes). Our analysis (see SI Appendix for details) indicates that the dominant contribution is the vdW interaction with the substrate $$ K_{eff}\approx \hspace{0.17em}{V_{vdW}}^{\text{'}\text{'}}\left(h_0\right)$$. The quantity “effective stiffness” does not have the conventional units of other measurables that are used to characterize stiffness in contact mechanics. Instead, we use the term “effective stiffness” as it is typically employed in the elasticity literature, where it refers to a “spring constant” associated with the resistive force exerted by a unit area of a compliant substrate. Hence the units are $$ \frac{\raisebox{1ex}{$\left[Force\right]$}\!\left/ \!\raisebox{-1ex}{$\left[Length\right]$}\right.}{\left[Area\right]}=\hspace{0.17em}\frac{\left[Energy\right]}{{\left[Length\right]}^4}$$. In our analysis, this term encapsulates the substrate resistance to the formation of ripples in the layer above it (thereby setting the wavelength through a balance with the bending rigidity of the sheet). In implementing this concept in our study, the relevant energy scale is not the minimum of the substrate-layer attraction energy/area, but rather the penalty for deviations from the minimal value (i.e., second derivative of $$ \frac{Energy}{Area}\to \frac{Energy}{Lentg{h}^4}$$).

Having established the nature of the effective stiffness in Eq. 2, let us now address the bending modulus B. Generally, there are two regimes for the bending of multilayer 2D crystals. At low bending curvature (large $$ \ell $$) there is no slippage between the planes, and the general continuous mechanics can be applied. For high bending curvature (low $$ \ell $$), when the layers are allowed to slip with respect to each other, the layers behave independently, so the bending rigidity is simply given by the sum of the individual bending rigidities of individual monolayers comprising the membrane. We define an effective bending rigidity, $$ B_{eff}$$, as:

where N is the number of layers, λ and μ are in-plane elastic constants (Lamé coefficients), d is the distance between layers, B is the bending rigidity of a single layer, and $$ {\text{Γ}}_s$$ determines the shear mode of multilayered samples. At long wavelengths, the mutual slippage between the layers is suppressed, and the effective bending rigidity scales as the cube of the number of layers, being determined by the in-plane bulk modulus, in agreement with the general theory of elastic slabs in three dimensions (33). At short wavelengths, intralayer slippage is prominent, such that the mechanical response of each layer is practically independent on other layers, and consequently, the bending rigidity equals to the number of layers multiplied by the bending rigidity of a single layer. The crossover wavelength $$ {\ell }^{\ast }$$ is determined by the shear stiffness and takes place at wavelengths and energies that belong to the lowest energy mode with out-of-plane modulations. The derivation of Eq. 3 is given in SI Appendix.

Fig. 4 D and E show the wavelength as a function of the number of layers for hBN and MoS₂, respectively, and fits of the experimental points to the “local wavelength law” (Eq. 2), with substrate-dominated stiffness. For these fits, the bending rigidity has been considered as the effective bending rigidity above the crossover length from Eq. 3; thus it scales as the cube of the number of layers, making $$ \ell $$ ∝ N^3/4. Note that Eq. 3 implies that above the crossover $$ {\ell }^{\ast }$$ the wavelengths tend to zero for n = 1, and hence the fits are valid for n > 1 only. We have used the Lamé coefficients λ and μ for hBN and MoS₂ given by refs. 343536–37 (see SI Appendix, Table S1 for a compilation of the employed elastic constants). From the fits we obtain the effective substrate stiffness for each material, resulting in $$ K_{eff}$$(hBN) = 2.8 × 10⁻⁷ eV Å⁻⁴ and $$ K_{eff}$$(MoS₂) = 5.1 × 10⁻⁷ eV Å⁻⁴. Using these values for the effective substrate stiffness we can obtain the expected wavelengths for the case n = 1, substituting the monolayer bending rigidities for each material into Eq. 2 (black diamonds in Fig. 4 D and E, compatible with the experimental observed wavelengths). Fig. 4F summarizes the results for both materials, where we plot the wavelength as a function of $$ B_{eff}^{1/4}$$ and fit the data to Eq. 2. We find excellent agreement between the theory and the experimental observations.

Let us turn now to commensurate configurations, where the hBN membrane is crystallographically aligned with the underlying graphene membrane. In this case, the behavior is remarkably different. Small rotation angles between the lattices of the two crystals form a periodic structure with different stacking configurations, yielding a moiré pattern (38). Fig. 5 corresponds to a monolayer of hBN aligned on top of a graphene flake presenting bubbles, as in the incommensurate states. It can be observed, however, that in the commensurate case the wrinkles exhibit a substantial deviation from radial orientation and become constrained to the moiré periodicity originated by the strong interaction between the top and the bottom layers. The wavelength of the wrinkles in this case is $$ \ell $$ ∼ 7.8 nm, the same as the moiré periodicity, and they follow preferential 60° directions, presenting abrupt kinks where it is necessary to accommodate the underlying moiré pattern.

This higher interaction is directly reflected in the aspect ratio of the bubbles. Whereas the typical H/R ratio for the incommensurate state of a monolayer hBN on graphene is 0.11 ± 0.02, for the commensurate state this value increases to H/R = 0.16 ± 0.02 (SI Appendix, Fig. S4). As mentioned in the introduction, the aspect ratio H/R is proportional to the ratio $$ {\left(\frac{\text{Γ}}{Y}\right)}^{1/4}$$, where Γ is the membrane-substrate surface energy and $$ Y\gg \text{Γ}$$ is the 2D membranal stretching modulus. In particular, from the measured values of the aspect ratio, H/R, and using a Young modulus, Y = 289 N m⁻¹ = 18 eV Å⁻², for monolayer hBN (3), we can obtain the vdW interaction, $$ \mathrm{\Gamma }$$, from $$ \frac{H}{R}=C_0{\left(\frac{\mathrm{\Gamma }}{Y}\right)}^{1/4}$$ (9), for both incommensurate and commensurate cases, resulting in $$ {\mathrm{\Gamma }}_{incomm}=$$ 0.003 eV Å⁻² and $$ {\mathrm{\Gamma }}_{comm}=$$ 0.013 eV Å⁻² (Fig. 3B). The value of $$ {\mathrm{\Gamma }}_{incomm}$$ is in line with other estimates of the vdW energy in similar setups (9). The membrane-substrate surface energy in the commensurate state is ∼4 times higher than for its counterpart in the incommensurate state (see Fig. 3B, which shows schematics of the vdW potential for the commensurate and incommensurate cases, highlighting the pronounced difference in their respective profiles). This stronger interaction is likely to drastically decrease the presence of water molecules (outside the bubble) between the membrane and the substrate, in comparison to the incommensurate state. This enhanced interaction is also reflected in the fact that the wrinkles in the commensurate state present almost atomic-scale heights (Fig. 5 B and C). Carrying out a similar analysis as for the incommensurate case, our experimental observations in the commensurate state indicate that the effective stiffness $$ K_{eff}$$, and, consequently, the vdW interaction with the substrate, is substantially enhanced in the commensurate case. From Eq. 2 we can calculate $$ K_{eff}$$ for this case, resulting in $$ K_{eff}$$(hBN)_comm ∼ 4 × 10⁻⁵ eV Å⁻⁴, which is two orders of magnitude higher than for the incommensurate state (please note the different steepness of the vdW potential in Fig. 3B). This may not be surprising, since the presence of water and other molecules in the incommensurate wrinkles may “soften” the interaction. The higher interaction in the commensurate case, reflected in a significantly higher magnitude of the vdW attraction and the effective stiffness, suggests that commensurate moiré structures of lengths on the order of 10 nm or larger are very stable.

Our experimental observations portray a nontrivial, multiscale interplay of the vdW substrate-membrane attraction, the high membranal resistance to stretching, and its low bending rigidity, driven by the presence of trapped substances between the membrane and the substrate. The primary response, which has been reported previously and is barely affected by bending rigidity, is the aggregation of trapped substances into highly pressurized bubbles of radius R, thereby causing delamination of the membrane from the substrate at r ≪ R (where h(r) > h_max) in Fig. 3, and, consequently, a large radial tension, $$ {\sigma }_{rr}\sim {\text{Γ}}^{1/2}{Y}^{1/2}\gg \text{Γ}$$ in the delaminated portion of the membrane. The primary response, which is evidently governed by the membranal stretching resistance (through the Young’s modulus Y) and the overall magnitude of the vdW attraction ($$ \text{Γ}\approx V\left(h_0\right)$$), underlies a secondary effect, due to the highly anisotropic stress in the solid membrane, whereby a hoop compression induced by the radial tension in the vicinity of the bubble is suppressed with the aid of radial wrinkles. The consequent suppression of the elastic energy is governed by much lower contributions, associated with the membranal bending modulus B and the steepness of the vdW potential, $$ K_{eff}\approx V\text{'}\text{'}\left(h_0\right)$$. Thus, measurements of the various features of the deformed membrane can be harnessed to provide a very convenient, indirect metrological probe. In the current study, we employed this system to determine the dependence of the bending rigidity on the number of atomic layers in the top membrane and to examine the strong effect of commensurability and humidity on the magnitude and steepness of the vdW potential. We expect that future studies will advance further this approach for metrology of vdW heterostructures.