Van der Waals interaction affects wrinkle formation in two-dimensional materials

Significance Instabilities and formation of complex patterns play an extremely important role in many branches of science. Understanding the origin of instabilities is very important for many technologies. We demonstrate the formation of instabilities (radial wrinkles around bubbles) in atomically thin crystals. Two-dimensional crystals, arranged into van der Waals heterostructures, allow full control over the mechanical properties of the components and the interaction between them, making such experiments easy to reproduce and allowing one to understand the basic, underlying physics of such effects. Such heterostructures allow the observation of instabilities around structural defects (bubbles) where the crystal structure plays a significant role, providing insights into these effects. Such instabilities allow assessment of the elastic and interaction properties of the membranes.

Height scales: 3 nm for 1L and 10 nm for the rest of the cases. The gray squared-dotted lines correspond to the substrate whereas the pink circled-dotted lines to wrinkles. The arrows in the profiles indicate when the bubble reaches the substrate level in the areas without wrinkles. In all the cases, the height of the wrinkles decays exponentially, such that most of the wrinkled zone is nearly flat in the radial direction and smaller than 1 nm in the majority of its length.

Substrate-induced stiffness of a membrane on a supporting substrate
The effective stiffness ( ) can originate from multiple contributions (5):  4 ) acts to suppress the defection (amplitude) of wrinkles from the rest (naturally planar) state of the membrane. Generally, the effective stiffness, ( ), comprises three types of restoring (amplitude-suppressing) forces: i. A substrate-induced stiffness, , attributed to an energetic cost, , emanates from the normal force exerted on the membrane by a supporting substrate, = (ℎ ). For polymer sheets, such a normal force reflects the energetic cost required to deform the substrate (e.g. Young's modulus of a compliant solid substrate (6), or the g.p.e. of a liquid bath). In our study, we can consider the supporting substrate as infinitely rigid and hence follow Zhang-Witten (7), assuming that the only normal force exerted by the substrate on the membrane originates from the change in the vdW energy. Namely, • = (ℎ ) • is the force required to bring a unit area of the membrane to a distance ℎ = ℎ + from the substrate (where ℎ = ℎ is the minimal-energy distance of a flat membrane from the substrate). Hence, = (ℎ ). ii. A tension-induced stiffness, ( ), is attributed to an energetic cost, ∝ ′( ) , of radial variation of the wrinkle amplitude, ( ), in the presence of radial tension ~Γ "/ $ "/ , along wrinkles (8). Since the radial variation of the wrinkle amplitude occurs over a scale ∝ %, we may estimate, , is attributed to an energetic cost of the strain associated with deforming a curved "envelope" shape, ' ( ), with radial curvature ~ ' ( ), due to the Gaussian curvature induced by superimposing on it azimuthal undulations (5).
Our experimental observations enable us to rule out the relevance of ( ) and ( ). First, observation (c) in the main text (the non-oscillatory component of the wrinkles decays exponentially, such that most of the wrinkled zone is nearly flat in the radial direction, and highly curved in the azimuthal direction) implies that, while ( ) may be significant very close to the bubble, it decays exponentially rapidly, and does not affect the membrane in the vast part of the wrinkled zone. Second, observation (e) (for a given type and number of layers of the top membrane, the average wavelength ℓ does not depend on the size (radius or height) of bubbles) implies that ( ) is also negligible, otherwise the average wrinkle wavelength would exhibit a strong dependence on the bubble's radius (ℓ ∝ √%). Finally, observation (b) (the wrinkle number +( ) increases with radial distance r) indicates a spatially-uniform wavelength, consistently with ≈ .
A straightforward relation ~ -. /0 1 could be used to estimate a characteristic width of the vdW potential well. This expression leads to length scales of the order of 1 nm for the commensurate case, and 10 nm for the incommensurate one. These estimates, which do not seem realistic, will be significantly reduced if one considers, as seems likely, that the very amplitude of the wrinkles and the thermal fluctuations of the membranes at room temperature smooth the vdW potential.

Bending rigidity of stacks of weakly coupled layers
The bending rigidity of isotropic, or almost isotropic, elastic slabs is proportional to the in-plane bulk modulus, and it increases as the cube of the thickness of the slabs (9). A stack of layers coupled by the van der Waals interaction has highly anisotropic elastic constants, and a different response to bending. The coupling between neighboring layers, at long wavelengths, can be described by three force constants, which determine i) the energy required for a change in the interlayer distance, that is, the breathing mode, ii) the energy required for a lateral displacement of one layer with respect the other, the shear mode, and iii) the coupling between a compression within the layer and a change in the interlayer distance, which determines the out of plane Poisson ratio.
The coupling between in plane deformations and bending in a slab arises from the shear deformation which arises when the slab is bent. Hence, it vanishes when the layers can freely slide one over the other. Neglecting the third coupling mentioned above, the elastic energy of a slab of 2 layers can be written as (10): 3 = ∑ 5 . ⃗ 7 8 9: ; < ; 0 + : = < = 0 > + ? @(: ; < ; 0 ) + 9: = < = 0 > + " 9: ; < = 0 + : = < ; 0 > AB + 0 Where V and ? are in plane elastic constants, Γ and Γ describes the interlayer breathing and shear stiffness, Z is the bending rigidity of each layer, and . is the distance between layers. The bending rigidity of a multilayer with N layers can be understood by analyzing the normal modes defined using the energy in eq. (2).
We focus on the low energy modes, described by smooth variations of the displacements between adjacent layers. A simple ansatz which is consistent with the continuum limit of isotropic slabs, which describes normal modes of eq. (2) is: with ^=^;, and ℎ and < are the amplitudes of out of plane and in plane displacements. The deformation described in eq. (3) does not require a change in the distance between layers, so that it does not depend on the out of plane elastic stiffness, Γ . The insertion of eq. (3) into eq. (2) leads to two modes for each value of ^, depending on the relative sign between ℎ and < . We consider the mode of lowest energy, e f (^), whose dispersion is: Where A is the area of the unit cell, g = 2 × 2 × g p is the mass of the unit cell, and g p is the mass of the carbon atom. Expanding eq. (4), we obtain: Which allows us to define an effective bending rigidity: where the wavelength is ℓ = (2u)⁄ . At long wavelengths, the effective bending rigidity scales as the cube of the number of layers, and it is determined by the in plane bulk modulus, in agreement with the general theory of elastic slabs (9). At short wavelengths, the bending rigidity scales like the number of layers, and it is proportional to the rigidity of a single layer. The crossover is determined by the shear stiffness.
We can use the paradigmatic case of graphene to gain some insight on the values of this crossover length. The elastic constants of single layer graphene are: And the interlayer distance is . ≈ 3.5 Å. The parameter Γ determines the shear mode of graphite and multilayered samples, where = 93√3• > 2 ⁄ is the area of the unit cell of graphene, and • ≈ 1.4 Å is the distance between neighboring carbon atoms.
The frequency of the shear mode in bilayer graphene has been calculated and measured in Ref. (11), and a related calculation for multiwalled nanotubes can be found in Ref. (12). The Raman frequency of the interlayer shear mode in multilayered samples obtained in Ref. (11) depends on the number of layers, 2. It ranges between e f ≈ • × 31 cm 1" for a bilayer and e f ≈ • × 44 cm 1" for graphite, where • is the velocity of light. The dependence of the shear frequency on the number of layers implies the existence of couplings between layers that are not nearest neighbors, not included in eq. (2). Neglecting this effect, we obtain the value of Γ , and a crossover length, ℓ * between the two regimes described in eq. (6): Γ ≈ 17 meV Å 1 ℓ * (2) = m k" l (b j 1b) l (8N s) " (b1")C K (9) This length increases as the bulk modulus of a monolayer increases, and it decreases as the interlayer shear stiffness increases. For thick multilayers, we find lim b→ˆℓ * (2) ≈ 240 × 2 nm , which is more than two orders of magnitude larger than the width of the stack, ‰ = (2 − 1) × .. For a graphene bilayer, the length is ℓ * ≈ 55 nm.
In very rigid monolayers, the sliding of the monolayers is favored, and the bending rigidity of the stack is determined by the bending rigidity of each layer. On the other hand, in isotropic systems Γ ~ V + 2?, and ℓ * ~2 × .. The bending rigidity is independent of the bending rigidity of the individual layers for length scales larger than the width of the slab.
It is interesting to note that in twisted bilayers with a small twist angle most atoms in different layers are not in registry, and the interlayer shear is very small or it even vanishes, when the resulting Moiré pattern is incommensurate with the atomic lattices (12). Then, the bending rigidity of a small angle twisted graphene bilayer is always twice the bending rigidity of each layer.