Tunable elastic wave bandgaps by strain engineering of multilayered van der Waals metamaterials

Abstract

Multilayered van der Waals (vdW) metamaterials exhibit exceptional mechanical and thermal properties, with the ability to tune these characteristics based on material composition, stacking sequences, and the number of layers. Due to the extremely small scale of these structures, the influence of vibration and noise becomes significant. This study uses molecular dynamics (MD) simulations to investigate the bandgap characteristics of elastic waves in multilayered h-BN/MoS₂ vdW metamaterials. It demonstrates that the position and width of the broadband terahertz elastic wave bandgap in these materials can be adjusted by modifying the longitudinal strain and stacking order. The bandgap is notably sensitive to the rate of longitudinal strain, with its frequency increasing under compressive strain and decreasing under tensile strain. The effects of tensile and compressive strains on the bandgap are asymmetric; compressive strain can notably alter the bandgap width and even cause it to disappear. Additionally, the bandgap can be manipulated by shifting the position of the h-BN atomic layer or by rotating the h-BN layer by 90° within the metamaterial's unit cell. This research offers a novel method for tuning elastic wave propagation in multilayered vdW heterostructures.

---

1. Introduction

Multilayered van der Waals (vdW) metamaterials, formed by vertically stacking two-dimensional (2D) materials through vdW forces, exhibit outstanding optical,^1–3 electrical,^4–6 mechanical,^7–10 and thermal^11,12 properties. These properties can be tuned by altering the material type, stacking order, twist angle, and applied strain. Typically, multilayered vdW metamaterials can be as thin as 10 nm and outperform the individual materials of which they are composed.^13–17 These advantages have spurred significant theoretical and experimental research into vdW heterostructures for optoelectronic devices.^18–24

The impact of vibrations in nanostructures becomes increasingly significant as devices shrink in size, and many studies have explored this topic.^25–30 However, research on the vibration behavior of multilayered vdW metamaterials has been limited due to experimental challenges. He et al.³¹ developed a continuum plate model to study the vibration of multilayered graphene sheets, based on a formula predicting vdW interactions between any two sheets. The vdW forces^32–34 are critical in influencing the vibrations of multilayered vdW structures. Zhang et al.³⁵ examined the vibration behavior of homogeneous multilayered structures using molecular dynamics (MD) simulations and a sandwich plate model. They found that the natural frequencies of multilayered black phosphorus (BP) and h-BN varied significantly with the number of layers, while changes in multilayered graphene were less pronounced. The sandwich plate model, which accounts for the shear modulus between adjacent layers due to vdW forces, more accurately predicted the vibrations of multilayered vdW metamaterials compared to MD simulations. Hou et al.³⁶ investigated the vibration characteristics of BP and graphene heterostructures using a mesh-free method. Their results indicated that, while vdW forces between arbitrary and adjacent layers made little difference in low-order natural frequencies, the difference became significant in higher-order frequencies.

In addition to the number of layers, several other factors can influence the vibration behavior of multilayered vdW metamaterials. Homogeneous multilayered vdW materials typically do not exhibit vibration bandgaps. However, few-layer graphene (FLG) with two highly tensioned outermost layers can develop an elastic wave bandgap in the terahertz range, resulting in significant vibration isolation.^37,38 Strain can significantly alter the mechanical properties of two-dimensional materials.^39,40 The impact of stress on elastic waves has also been observed in two-dimensional structures. Kirchhof et al.⁴¹ found that applying out-of-plane electrostatic gates to monolayer suspended graphene phonon crystals caused the quasi-bandgap to shift upward with increasing pressure, along with a corresponding increase in sound velocity. As biaxial pressure (σ_xx/σ_yy = 1) increases, the bandgap rises and its range expands. In contrast, with increasing uniaxial tension (σ_xx/σ_yy > 1), the bandgap frequency shifts upward, but its range gradually decreases until it disappears.⁴²

One of the most effective strategies for engineering interlayer interactions is through stacking manipulation, particularly by adjusting the twist angle^16,43–47 and applying mechanical strain.^48–50 Terahertz phonon engineering plays a crucial role in phonon–optical coupling and high-frequency thermal management. Yoon et al.⁵¹ have demonstrated the generation of broadband terahertz phonons up to 3 THz using FLG, the sensitive detection of terahertz phonons with monolayer WSe₂, and the flexible manipulation of terahertz waves through carefully designed h-BN/WSe₂/FLG vdW heterostructure stacks. They also modeled the observed pump–probe signals and corresponding atomic motions using a one-dimensional mass-spring chain model. Additionally, our previous work⁵² introduced a mass-spring model-based approach to describe terahertz elastic wave propagation in multilayered h-BN/MoS₂ vdW metamaterials, investigating how different 2D materials and stacking sequences influence the bandgap frequency. This method is versatile and applicable to other 2D layered materials as well. In-plane uniaxial tensile strain has a minimal impact on phonon transport perpendicular to the cross-sectional direction in multilayer heterojunctions.⁵³ However, applying tensile strain along the cross-section direction significantly suppresses out-of-plane atomic vibrations, thereby hindering phonon transport across the heterojunction, while compressive strain enhances phonon transport. Nevertheless, excessive strain leading to structural folding can substantially reduce phonon transport. Mechanical strain can arise from lattice mismatches in multilayered vdW structures and the application of external loads. To our knowledge, no prior research has explored the regulation of elastic wave bandgaps in multilayered vdW metamaterials through stacking orders and strain applied along the stacking direction.

The MD method is an effective approach for studying the vibration of nanostructures,^44,54–57 especially when the experimental operations are difficult, and an effective theoretical method is lacking. This work examines the bandgap characteristics of elastic waves in h-BN/MoS₂ vdW metamaterials using MD simulations. This paper is organized as follows: Section 2 details the molecular model of the multilayered h-BN/MoS₂ vdW metamaterial and the MD simulations used to determine its bandgap characteristics. The effects of longitudinal strain and stacking orders on bandgap characteristics are analyzed in Sections 3 and 4, respectively. Finally, concluding remarks are presented in Section 5.

2. Molecular dynamics model and simulations

MD simulations are utilized to investigate the bandgap characteristics of multilayered h-BN/MoS₂ vdW metamaterials. The molecular model of the metamaterial is first constructed, with a perspective view illustrated in Fig. 1(a). The unit cell of the metamaterial comprises vertically stacked single layers of h-BN and MoS₂ atoms, and the metamaterial itself consists of 5 such unit cells. The lattice constants for h-BN and MoS₂ are a_B = 2.51 Å and a_M = 3.185 Å, respectively. Both the h-BN and MoS₂ nanosheets measure 5.9 × 5.9 nm², and the entire system includes 12 990 atoms. The front and side views of the unit cell, as depicted in Fig. 1(b), provide a detailed representation of the stacking order.

Fig. 1 Molecular model of the multilayered h-BN/MoS₂ vdW metamaterial. (a) Perspective view of the metamaterial. (b) Front view and side view of the unit cell of the metamaterial. (c) Boundary conditions of the metamaterial, displacement excitation is applied to the lowest atomic layer, and the uppermost atoms are fixed. The red dot near the uppermost of the metamaterial is the calculation point of the transmissibility.

All MD simulations are conducted using a large-scale atomic/molecular massively parallel simulator (LAMMPS).⁵⁸ The bonding interactions within the Mo–S system are described using the reactive empirical bond-order (REBO) potential parameters developed by Liang et al.^59,60 For the in-plane bond interactions within the h-BN layers, the Tersoff potential⁶¹ is employed. The interlayer interactions between the different layers are modeled using the Lennard-Jones (LJ) potential. The energy of the LJ potential is expressed as follows:

where r is the distance between interacting atoms, and ε and σ are the well depth and the equilibrium distance, respectively. The LJ potential parameters for Mo–Mo, and S–S are taken from ref. 62 and B–B, N–N and B–N are taken from ref. 63, as shown in Table 1. The remaining interatom interaction parameters shown in Table 1 are obtained using the Lorentz–Berthelot mixing rules:A time step of 1 femtosecond (fs) is employed, with periodic boundary conditions applied in all directions. To simulate a finite thickness model, vacuum layers of 20 Å are added at both ends along the z-direction, which corresponds to the stacking direction of the 2D materials. This distance is sufficient to prevent interactions between periodic images, given a cutoff radius of 10 Å. The system undergoes energy minimization using the steepest descent (SD) algorithm and the damped dynamics method of fast inertial relaxation engine (FIRE).⁶⁴ This process is repeated until the system reaches its minimum potential energy. Subsequently, simulations are carried out in the microcanonical ensemble. Sinusoidal displacement excitations with a magnitude of A_in = 0.005 Å at various frequencies are applied to the lowest atomic layer in the direction perpendicular to the layers, and the uppermost atoms are fixed during MD simulations, as depicted in Fig. 1(c).

Table 1 The LJ potential parameters of h-BN and MoS₂

	S–S	Mo–Mo	B–B	N–N	B–N	Mo–S	B–Mo	B–S	N–Mo	N–S
ε (meV)	13.86^62	0.58595^62	4.116^63	6.281^63	5.085	2.8498	1.553	7.553	1.918	9.330
σ (Å)	3.13^62	4.2^62	3.453^63	3.365^63	3.409	3.665	3.8265	3.2915	3.7825	3.2475

---

The sinusoidal displacement excitation with frequency ω is applied to the metamaterial, and then the time-domain signal x_n of the oscillation of the red dot in Fig. 1(c) is collected in the range of t_total = 0–1 ns. x_n represents the value of the signal at t = nΔt(n = 0, 1, 2,…, N–1), where N is the number of sample points, and Δt = t_total/N is the sampling interval. The signal in the frequency domain can be obtained by applying a fast Fourier transform (FFT) to the time domain signal. The power P_ω, that is, the power corresponding to the frequency component ω in the total power of atomic oscillations near the end of the metamaterial as shown by the red dot in Fig. 1(c), is equal to the square of the amplitude of the frequency component ω in the frequency domain signal. Finally, the transmissibility at this excitation frequency ω can be obtained by T = 10 log₁₀(P_ω/P_in), where P_in is the input power. By changing the excitation frequency ω and repeating the above calculation process, the transmissibility curve of the metamaterial as shown in Fig. 2(a) can be obtained. The continuous frequency region of T < −27 dB is defined as the bandgap region.

Fig. 2 Bandgap characteristics of the multilayered h-BN/MoS₂ vdW metamaterial. (a) Transmissibility of the metamaterial with fixed and free boundary conditions. The blue dots A and B at 1 THz and 2 THz, located inside and outside the bandgap region, respectively, are the specified frequency points to be further studied later. (b) Displacement time history of the metamaterial from 0.4 ns to 0.5 ns when the excitation frequencies at points A and B, respectively. The black and red curves are the displacement of atomic oscillations of excitation (x_in) and the response of the red dot in Fig. 1(c) (x_n), respectively. The inset shows the displacement time history of the metamaterial from 0.45 ns to 0.46 ns at 2 THz.

The red line in Fig. 2(a) clearly indicates that the vdW metamaterial with a fixed boundary exhibits an elastic wave bandgap in the terahertz range. Frequency points A and B at 1 and 2 THz, corresponding to the blue dots outside and inside the bandgap in Fig. 2(a), are selected for further analysis. The displacement time histories of the vdW metamaterial from 0.4 to 0.5 ns under displacement excitations at these two points A and B are shown in Fig. 2(b). The black curve represents the excitation displacement (x_in, x_in is the time-domain signal of the displacement excitation), and the red curve indicates the response displacement (x_n) of the red dot in Fig. 1(c). The results show that elastic waves within the bandgap range are significantly suppressed. The effect of damping is not considered in the MD simulations, which leads to the existence of amplitude components of other frequencies. Therefore, the value of the transmissibility T at 2 THz is inconsistent with the attenuation effect reflected by the time domain signal.

In addition, transmissibility T of the vdW metamaterial with unfixed uppermost atoms under sinusoidal displacement excitation is also calculated, and the results are shown in Fig. 2(a) by the black dashed line. There is no significant difference between the bandgap ranges obtained when the uppermost atoms of the metamaterial are fixed and unfixed during MD simulations. Following this, the effects of longitudinal strain and stacking order on the bandgap characteristics are investigated.

3. Effects of longitudinal strain on the bandgap characteristics

To investigate the effect of longitudinal strain on the bandgap characteristics of multilayered vdW metamaterials, simulations were conducted to observe the propagation of elastic waves under both tensile and compressive strains applied along the stacking direction. The strain is expressed as ε = (l – l₀)/l₀ × 100%, where l₀ is the original length of the metamaterial along the stacking direction, and l is the length after strain is applied. A positive ε indicates tensile strain, while a negative ε indicates compressive strain. In the simulation, the uppermost atomic layer is fixed, and the lowermost atomic layer is shifted incrementally, minimizing the system's energy with each shift until the desired ε value is reached. Finally, the transmissibility curve for a specific ε value is obtained by fixing the uppermost atomic layer and applying displacement excitation with varying frequencies to the lowermost atomic layer.

Fig. 3 illustrates how the bandgap is influenced by both compressive and tensile strain. The width and frequency of the metamaterial's bandgap are closely related to the strain (ε). Initially, as the strain increases, the bandgap widens. However, with further strain, the bandgap gradually narrows. Notably, the effects of compressive and tensile strains on the bandgap are asymmetric. The bandgap changes more rapidly under compressive strain compared to tensile strain. Under compressive strain, the bandgap quickly diminishes and eventually vanishes when the strain ε is −4.2%. Conversely, under tensile strain, the bandgap narrows at a much slower rate, and it disappears only when the strain reaches a certain threshold. As compressive strain ε increases beyond 4.2%, the critical frequencies between the passband and bandgap continue to rise.

Fig. 3 Bandgap changing induced by compressive and tensile strains. Red and blue dots represent the upper and lower boundaries of the bandgap, respectively. The inset is the curve of the interlayer energy of the unit cell of the h-BN/MoS₂ vdW metamaterial changing with the spacing between layers.

The inset in Fig. 3 depicts the energy-spacing curve for the h-BN and MoS₂ atomic layers. As strain increases, the layer spacing increases, leading to an increase in the vdW force. However, the increase in the vdW force is more pronounced under compressive strain compared to tensile strain. Consequently, as compressive strain increases, the vdW force rises sharply, causing significant changes in the bandgap position and frequency. This results in elastic waves traveling more easily between the atomic layers, shifting the bandgap frequency upwards and narrowing its width. In contrast, as tensile strain increases, the layer spacing increases, making it more difficult for elastic waves to propagate between the atomic layers. This causes the bandgap position to shift slowly downwards.

Fig. 4(a)–(d) present the transmissibility of the metamaterial under strain ε of 2.4%, 7.51%, −2.4%, and −4.2%, respectively. The transmissibility exhibit significant variation depending on the strain applied. Fig. 4(d) reveals that the disappearance of the bandgap does not imply that elastic waves can propagate through the metamaterial at all frequencies; instead, it indicates that the bandgap is immediately followed by the passband.

Fig. 4 Transmissibility of the vdW metamaterial at different strains. (a)–(d) The transmissibility of the metamaterial at strains of 2.4%, 7.51%, −2.4% and −4.2%, respectively. Points with different colors in the figure corresponding to points in Fig. 3.

To investigate the vibration behavior of the metamaterial at various strain and frequencies, displacement time histories corresponding to the different strains and frequencies marked by colored dots in Fig. 4 are presented in Fig. 5(a)–(d). The results show that elastic waves within the bandgap frequency are significantly suppressed. As depicted in Fig. 5(b), when ε = 7.51% and f = 0.1 THz, the vibration transitions into a new state of steady vibration at a certain point. When the amplitude of the oscillation is too large, the interaction between the first and second atomic layers exists but is greatly weakened, and the spacing between the two layers suddenly increases significantly. Since the ends of the vdW metamaterials are fixed, the other atomic layers do not completely separate from the first layer. The upper half of the metamaterial vibrates with a large amplitude as a whole, while the first atomic layer still vibrates around the initial position under displacement excitation. The detailed vibration process and displacement contours are shown in Fig. 6. Notably, the vibration only diverges at low frequencies, the amplitude becomes larger, and then the steady-state transition occurs. As the frequency increases, a range of frequencies where vibration reduction occurs also emerges. Furthermore, the vibration becomes unstable more rapidly as the strain increases, as illustrated in Fig. 5(e).

Fig. 5 Displacement time histories of the vdW metamaterial. (a)–(d) Displacement time histories of the metamaterial from 0 to 0.4 ns at different strains and frequencies, respectively. The green, blue, purple, orange and grey curves represent the displacement of corresponding colored dots in Fig. 4 and the excitation displacement. (e) Steady-state transition process with the increase of strain. Enlarged views correspond to the displacement time history within the boxes.

Fig. 6 Vibration of the vdW metamaterial with ε = 7.51% and f = 0.1 THz. (a) Partial structure of the metamaterial. (b) and (c) Displacement time histories of the oscillation of atoms in (a), and the right side are the enlarged views of the red boxes. (d)–(g) Displacement contours of the metamaterial at different times in the enlarged view.

As the strain continues to increase beyond approximately 7.8%, the spacing between the first layer (h-BN) and the second layer (MoS₂) suddenly expands. Concurrently, the vdW force between these layers diminishes, and the interlayer interaction gradually approaches zero as the strain increases. The detailed stretching process is presented in the Appendix.

Overall, adjusting the longitudinal strain in multilayered vdW metamaterials enables dynamic control over the frequency range and width of the bandgap, allowing for precise regulation of the bandgap size. This process is purely mechanical, requiring the application of appropriate strain without involving any chemical alterations. Furthermore, it preserves the original structure and material composition of the metamaterial, making it a straightforward and effective method for tuning its properties.

4. Effect of stacking orders on the bandgap characteristics

The stacking order in 2D materials like graphene, BP, and others has a significant impact on their vibrational properties^44,65 and electronic properties.^66,67 To understand the effect of stacking order on the elastic wave bandgap of multilayered vdW metamaterials, new metamaterials B and C were designed, as shown in Fig. 7(a). Metamaterial B was created by shifting the h-BN atomic layer along the y-direction by a_B/4 within the unit cell of metamaterial A, as described in Section 2, while metamaterial C was formed by rotating the h-BN atomic layer 90° clockwise. The front, side, and top views of the unit cells for metamaterials A, B, and C are presented in Fig. 7(a).

Fig. 7 Molecular model and interlayer interaction of metamaterials A, B and C. (a) Front, side and top views of the molecular model of the unit cell of metamaterials A, B and C, respectively. (b) Total interaction energy of unit cells of metamaterials A, B and C vary with layer spacing. (c)–(e) Interaction energy map of all atoms in the h-BN atomic layer with atoms in other layers and the z-direction displacement of atoms in the h-BN layer at 2 THz of metamaterials A, B and C, respectively.

In multilayered vdW metamaterials, the factors influencing the elastic wave bandgap include interlayer interaction forces, density, and the strain on the atomic layers. The interaction between layers in these materials is primarily governed by vdW forces.^31,68 To explore the impact of different stacking orders on the system's energy, an analysis was conducted considering only the vdW forces within the unit cell. The equilibrium interlayer distance in the h-BN/MoS₂ unit cell was determined by minimizing the system's energy. Following this, the h-BN atomic layer was fixed while the MoS₂ atomic layer was gradually displaced along the z-direction. After each displacement, the system's energy, considering only the interlayer interaction, was calculated. Fig. 7(b) illustrates how the system's energy changes with varying layer spacing. The equilibrium distances—where the energy is lowest—between the h-BN and MoS₂ atomic layers for metamaterials A, B, and C show minimal variation. Additionally, the interlayer energy across the three stacking orders exhibits no significant differences.

To provide a more detailed view of the interaction energy between atomic layers, the interaction energy of each atom in the h-BN layer with all atoms in other layers within the cutoff radius was calculated, excluding the in-plane interaction energy. Fig. 7(c)–(e) present the interaction energy and displacement vector for metamaterials A, B, and C, respectively. The images show, from top to bottom, the interlayer interaction energy map of atoms in the h-BN layer and the z-direction displacement of atoms in the h-BN layer at 2 THz. The significant inherent lattice mismatch ε_m = 23.7% between the MoS₂ and h-BN layers leads to the formation of moiré patterns following geometry optimization, as depicted in Fig. 7(a). Here, ε_m = 2(a_M − a_B)/(a_M + a_B). Metamaterials A and B display moiré patterns with a period of a_P1 = 1.24 nm, while metamaterial C exhibits moiré patterns with a period of about a_P2 = 0.235 nm. Compared to metamaterial A, the energy distribution in metamaterial B shows a shift of a_P1/4 in the positions of the energy valleys and peaks along the y-direction, which corresponds to the shift of the h-BN atom layer by a_B/4. However, the energy distribution in metamaterial C differs significantly from those in metamaterials A and B. In metamaterial C, the energy peaks are more densely packed than in metamaterials A and B. The interaction energy corresponds closely with the moiré pattern, and the positions of atoms with the largest displacement generally align with the peaks of the potential energy profile.

The distribution of atoms at the positions corresponding to the peaks and valleys in the potential energy profile is generally consistent across different stacking orders. As a result, the vdW energy curves for the three stacking orders shown in Fig. 7(b) nearly overlap. However, while metamaterials A and B exhibit similar potential energy profile patterns, metamaterial C displays a distinct pattern. In metamaterial C, the maximum displacement is larger and more densely distributed compared to metamaterials A and B. Although the total interlayer interaction energy remains nearly identical across the stacking orders, the energy distribution and z-direction displacement differ significantly. Finally, the bandgap characteristics of metamaterials A, B, and C were examined using MD simulations, with the resulting transmissibility presented in Fig. 8. Despite the nearly identical interlayer energy among metamaterials A, B, and C, their bandgap characteristics differ significantly. In the bandgap range of metamaterial A, the atomic displacement in metamaterials A and B is relatively small, whereas in metamaterial C, the displacement is much larger, as illustrated in Fig. 7(c)–(e). The bandgap width increases when the h-BN layer is translated (metamaterial B), while the bandgap disappears entirely when the h-BN layer is rotated 90° (metamaterial C). Additionally, the bandgap characteristics of metamaterials D and E, which are formed by translating the h-BN layer by distances of a_B/3 and a_B/2 along the y-direction within the unit cell, were also investigated. Details of these findings are shown in Fig. 9.

Fig. 8 Transmissibility of metamaterials A, B and C. The different colored lines represent the transmissibility of metamaterials A, B and C as shown in Fig. 7(a), respectively.

Fig. 9 Energy maps and transmissibility of metamaterials A, D and E. (a) Molecular models and interlayer interaction energy maps of metamaterials A, D and E, respectively. (b) Transmissibility of metamaterials A, D and E. The different colored lines represent the transmissibility of the metamaterial A, D and E, respectively.

In summary, translating the h-BN layer along the y-direction alters the bandgap region, while rotating the h-BN layer by 90° changes the moiré pattern's period, leading to the disappearance of the elastic wave bandgap. These changes in the bandgap are primarily due to the variations in the intralayer strain resulting from the altered stacking order.

5. Conclusions

In summary, the bandgap characteristics of elastic waves in multilayered h-BN/MoS₂ vdW metamaterials are investigated using MD simulations. These metamaterials exhibit an elastic wave bandgap in the terahertz range, which can be influenced by longitudinal strain and stacking order. Longitudinal strain significantly affects the bandgap frequency range, with the bandgap width initially increasing and then decreasing under both compressive and tensile strains. Compressive strain has a more pronounced impact on the bandgap, eventually causing it to disappear. The bandgap can also be altered by translating the h-BN atomic layer along the y-direction, and it vanishes when the h-BN layer within the unit cell is rotated by 90°. The changes in the bandgap due to intralayer strain are linked to the interlayer potential energy profiles of different stacking orders. This work provides valuable insights for vibration reduction and the regulation of elastic wave bandgaps in multilayered vdW heterostructure devices.

Data availability

The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of interest

There are no conflicts to declare.

Appendix

This appendix provides a detailed discussion of the longitudinal tensile strain process in multilayered vdW metamaterials, as simulated using the MD method. Initially, the energy of the entire molecular model is minimized. The uppermost atomic layer of the metamaterial is held fixed, while the bottommost atomic layer is incrementally stretched by 0.06 Å at each step, followed by another round of energy minimization. This process of stretching and energy minimization is repeated, allowing for the calculation of the system's strain and stress in the tensile direction. Periodic boundary conditions are applied within the plane of the atomic layers, with a time step set at 1 ps.

Fig. 10 presents the stress–strain curve of the vdW metamaterial, with inset images showing the displacement vectors at various strain stages. Due to the relatively weaker interlayer interactions compared to intralayer interactions, the maximum longitudinal stress during the tensile process reaches only 1.61 GPa. At the initial stage of tensile deformation, where the deformation is small, the stress–strain relationship is nearly linear, indicating that interlayer vdW forces can be approximated as linear.⁵² As the strain increases, the rate of stress increase slows down. When the strain ε reaches 7.8%, continued stretching of the first layer causes its attractive force on the second layer to become weaker than the forces exerted by other layers on the second layer. Consequently, during the energy minimization process, all layers except the first layer retract towards their original positions, resulting in a final spacing between the first and second layers of approximately 6.7 Å, as depicted in the displacement vector in Fig. 10. The interaction between the first and second layers becomes minimal. As strain continues to increase, stress changes occur more gradually, resembling the trend seen in the energy-spacing curve in the inset of Fig. 7(b). Eventually, the spacing between the first and second layers approaches a cutoff radius of 10 Å.

Fig. 10 The stretching process of the multilayered vdW metamaterial. The inset figures are the displacement vector of the metamaterial at the corresponding strain values.

Acknowledgements

This work is supported in part by the National Science Fund for Distinguished Young Scholars under Grants No. 11925205 and in part by the National Natural Science Foundation of China under Grant No. 51921003 and U2341230.

References

1. X. Hong, J. Kim and S.-F. Shi, et al., Ultrafast charge transfer in atomically thin MoS2/WS2 heterostructures, Nat. Nanotechnol., 2014, 9(9), 682–686.
2. P. Rivera, J. R. Schaibley and A. M. Jones, et al., Observation of long-lived interlayer excitons in monolayer MoSe2–WSe2 heterostructures, Nat. Commun., 2015, 6(1), 6242.
3. Y. Liu, N. O. Weiss and X. Duan, et al., van der Waals heterostructures and devices, Nat. Rev. Mater., 2016, 1(9), 16042.
4. X. Cui, G.-H. Lee and Y. D. Kim, et al., Multi-terminal transport measurements of MoS2 using a van der Waals heterostructure device platform, Nat. Nanotechnol., 2015, 10(6), 534–540.
5. Y. Cao, V. Fatemi and A. Demir, et al., Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature, 2018, 556(7699), 80–84.
6. Y. Liu, J. Guo and E. B. Zhu, et al., Approaching the Schottky-Mott limit in van der Waals metal-semiconductor junctions, Nature, 2018, 557(7707), 696–700.
7. Z. C. Huang, Z. Z. He and Y. B. Zhu, et al., A general theory for the bending of multilayer van der Waals materials, J. Mech. Phys. Solids, 2023, 171, 105144.
8. X. Qian, F. Li and H. Tian, Stacking effects on dynamic mechanical behavior of bilayer hexagonal boron nitride under impact, Thin Wall Struct., 2024, 200, 111866.
9. Y. Yang, J. Ma and J. Yang, et al., Graphene/h-BN hybrid van der Waals structures with high strength and flexibility: A nanoindentation investigation, Thin Wall Struct., 2024, 195, 111341.
10. W. G. Ouyang, O. Hod and M. Urbakh, Parity-Dependent Moire Superlattices in Graphene/h-BN Heterostructures: A Route to Mechanomutable Metamaterials, Phys. Rev. Lett., 2021, 126(21), 216101.
11. Y. Gao and B. Xu, van der Waals Graphene Kirigami Heterostructure for Strain-Controlled Thermal Transparency, ACS Nano, 2018, 12(11), 11254–11262.
12. K. J. Tielrooij, N. C. H. Hesp and A. Principi, et al., Out-of-plane heat transfer in van der Waals stacks through electron-hyperbolic phonon coupling, Nat. Nanotechnol., 2018, 13(1), 41–46.
13. C. R. Dean, A. F. Young and I. Meric, et al., Boron nitride substrates for high-quality graphene electronics, Nat. Nanotechnol., 2010, 5(10), 722–726.
14. X. D. Duan, C. Wang and J. C. Shaw, et al., Lateral epitaxial growth of two-dimensional layered semiconductor heterojunctions, Nat. Nanotechnol., 2014, 9(12), 1024–1030.
15. K. S. Novoselov, A. Mishchenko and A. Carvalho, et al., 2D materials and van der Waals heterostructures, Science, 2016, 353(6298), aac9439.
16. R. M. Liu, J. Y. He and J. C. Zhang, et al., Moire Tuning of the Dynamic Behavior of a Twisted Bilayer van der Waals Material Resonator, J. Appl. Mech., 2022, 89(12), 121001.
17. A. Castellanos-Gomez, X. Duan and Z. Fei, et al., van der Waals heterostructures, Nat. Rev. Methods Primers, 2022, 2(1), 58.
18. G. Giovannetti, P. A. Khomyakov and G. Brocks, et al., Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations, Phys. Rev. B, 2007, 76(7), 073103.
19. R. Q. Cheng, F. Wang and L. Yin, et al., High-performance, multifunctional devices based on asymmetric van der Waals heterostructures, Nat. Electron., 2018, 1(6), 356–361.
20. S. J. Liang, B. Cheng and X. Y. Cui, et al., van der Waals Heterostructures for High-Performance Device Applications: Challenges and Opportunities, Adv. Mater., 2020, 32(27), 1903800.
21. H. Wang, H. Guo and R. Guzman, et al., Ultrafast non-volatile floating-gate memory based on all-2D materials, Adv. Mater., 2024, 2311652.
22. Z. Zeng, Z. Tian and Y. Wang, et al., Dual polarization-enabled ultrafast bulk photovoltaic response in van der Waals heterostructures, Nat. Commun., 2024, 15(1), 5355.
23. M. Jiang, T. Zheng and J. Zhang, et al., Gate-Modulated Polarity Transition and Polarization-Sensitive Photodetection Enabled by Sandwiching Anisotropic GeSe in vdW Heterojunction, Adv. Opt. Mater., 2024, 12(16), 2303217.
24. T. Zhang, X. Guo and P. Wang, et al., High performance artificial visual perception and recognition with a plasmon-enhanced 2D material neural network, Nat. Commun., 2024, 15(1), 2471.
25. R. Liu and L. Wang, Stochastically driven vibrations of single-layered graphene sheets, Sci. China: Phys., Mech. Astron., 2012, 55(6), 1103–1110.
26. J. Jiang, L. Wang and Y. Zhang, Vibration of single-walled carbon nanotubes with elastic boundary conditions, Int. J. Mech. Sci., 2017, 122, 156–166.
27. Z. P. Xu and Q. S. Zheng, Micro- and nano-mechanics in China: A brief review of recent progress and perspectives, Sci. China:Phys., Mech. Astron., 2018, 61(7), 074601.
28. K. M. Liew, Z. Pan and L.-W. Zhang, The recent progress of functionally graded CNT reinforced composites and structures, Sci. China:Phys., Mech. Astron., 2019, 63(3), 234601.
29. J. Zhang, L. Wang and G. Tai, et al., Splitting of vibration mode in an imperfect submicron circular plate, Acta Mech., 2021, 232(4), 1729–1739.
30. Y. Jing and L. Wang, Elastic wave propagation in a single-layered hexagonal boron nitride metamaterial, J. Appl. Phys., 2022, 131(18), 185104.
31. S. Kitipornchai, X. Q. He and K. M. Liew, Continuum model for the vibration of multilayered graphene sheets, Phys. Rev. B, 2005, 72(7), 075443.
32. X. Q. He, S. Kitipornchai and K. M. Liew, Resonance analysis of multi-layered graphene sheets used as nanoscale resonators, Nanotechnology, 2005, 16(10), 2086–2091.
33. M. Aydogdu, Vibration of multi-walled carbon nanotubes by generalized shear deformation theory, Int. J. Mech. Sci., 2008, 50(4), 837–844.
34. K. Kiani and H. Pakdaman, Nonlocal vibrations and potential instability of monolayers from double-walled carbon nanotubes subjected to temperature gradients, Int. J. Mech. Sci., 2018, 144, 576–599.
35. Y. Q. Zhang and L. F. Wang, Effects of the van der Waals Force on the Vibration of Typical Multi-layered Two-dimensional Nanostructures, Sci. Rep., 2020, 10(1), 644.
36. D. C. Hou, L. F. Wang and Y. Q. Zhang, Effects of van der Waals Forces on the Vibration of Stacked Multilayered Graphene/Black Phosphorus Heterostructures, Int. J. Struct. Stab. Dyn., 2021, 21(08), 2150115.
37. L. Lu, C. Ru and X. M. Guo, Metamaterial Vibration of Tensioned Circular Few-Layer Graphene Sheets, J. Appl. Mech., 2020, 87, 061009.
38. L. Lu, C. Q. Ru and X. M. Guo, Vibration isolation of few-layer graphene sheets, Int. J. Solids Struct., 2020, 185, 78–88.
39. M. A. Bissett, M. Tsuji and H. Ago, Strain engineering the properties of graphene and other two-dimensional crystals, Phys. Chem. Chem. Phys., 2014, 16(23), 11124–11138.
40. G. Barik and S. Pal, Strain-engineered BlueP–MoS2 van der Waals heterostructure with improved lithiation/sodiation for LIBs and SIBs, Phys. Chem. Chem. Phys., 2020, 22(3), 1701–1714.
41. J. N. Kirchhof, K. Weinel and S. Heeg, et al., Tunable graphene phononic crystal, Nano Lett., 2021, 21(5), 2174–2182.
42. J. N. Kirchhof and K. I. Bolotin, Mechanically-tunable bandgap closing in 2D graphene phononic crystals, npj 2D Mater. Appl., 2023, 7(1), 10.
43. M. Andelkovic, S. P. Milovanovic and L. Covaci, et al., Double Moire with a Twist: Supermoire in Encapsulated Graphene, Nano Lett., 2020, 20(2), 979–988.
44. J. C. Zhang, R. M. Liu and L. F. Wang, Negative interlayer shear effect on a double-layered van der Waals material resonator, Phys. Rev. B, 2021, 104(8), 085437.
45. D. W. Kang, Z. W. Zuo and Z. W. Wang, et al., Bandgap engineering of stacked two-dimensional polyaniline by twist angle, Appl. Phys. Lett., 2021, 119(6), 061602.
46. A. Rothstein, C. Schattauer and R. J. Dolleman, et al., Band gap formation in commensurate twisted bilayer graphene/hBN moire lattices, Phys. Rev. B, 2024, 109(15), 155139.
47. J. Wang, F. Cheng and Y. Sun, et al., Stacking engineering in layered homostructures: transitioning from 2D to 3D architectures, Phys. Chem. Chem. Phys., 2024, 26(10), 7988–8012.
48. Z. H. Ni, T. Yu and Y. H. Lu, et al., Uniaxial Strain on Graphene: Raman Spectroscopy Study and Band-Gap Opening, ACS Nano, 2008, 2(11), 2301–2305.
49. S. Shabani, D. Halbertal and W. J. Wu, et al., Deep moire potentials in twisted transition metal dichalcogenide bilayers, Nat. Phys., 2021, 17(6), 720–725.
50. E. Barré, M. Dandu and S. Kundu, et al., Engineering interlayer hybridization in van der Waals bilayers, Nat. Rev. Mater., 2024, 9(7), 499–508.
51. Y. Yoon, Z. Lu and C. Uzundal, et al., Terahertz phonon engineering with van der Waals heterostructures, Nature, 2024, 631(8022), 771–776.
52. Y. B. Jing, L. F. Wang and Y. Q. Gao, Mass-spring model for elastic wave propagation in multilayered van der Waals metamaterials, Appl. Math. Mech.-Engl. Ed., 2024, 45(7), 1107–1118.
53. X. Wu and Q. Han, Phonon Thermal Transport across Multilayer Graphene/Hexagonal Boron Nitride van der Waals Heterostructures, ACS Appl. Mater. Interfaces, 2021, 13(27), 32564–32578.
54. R. Liu and L. Wang, Nonlinear forced vibration of bilayer van der Waals materials drum resonator, J. Appl. Phys., 2020, 128(14), 145105.
55. L. Wang and H. Hu, Flexural wave propagation in single-walled carbon nanotubes, Phys. Rev. B, 2005, 71(19), 195412.
56. L. Wang, H. Hu and R. Liu, Thermal Vibration of Carbon Nanostructures, in Handbook of Mechanics of Materials, ed. S. Schmauder, C.-S. Chen, K. K. Chawla, N. Chawla, W. Chen and Y. Kagawa, Springer, Singapore, 2019, pp. 421–481.
57. Z. Wu, R. Liu and N. Wei, et al., Unexpected reduction in thermal conductivity observed in graphene/h-BN heterostructures, Phys. Chem. Chem. Phys., 2024, 26(5), 3823–3831.
58. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys., 1995, 117(1), 1–19.
59. T. Liang, S. R. Phillpot and S. B. Sinnott, Parametrization of a reactive many-body potential for Mo-S systems, Phys. Rev. B, 2009, 79(24), 245110.
60. T. Liang, S. R. Phillpot and S. B. Sinnott, Erratum: Parametrization of a reactive many-body potential for Mo-S systems (vol 79, 245110, 2009), Phys. Rev. B, 2012, 85(19), 199903.
61. A. Kinaci, J. B. Haskins and C. Sevik, et al., Thermal conductivity of BN-C nanostructures, Phys. Rev. B, 2012, 86(11), 115410.
62. J. A. Stewart and D. E. Spearot, Atomistic simulations of nanoindentation on the basal plane of crystalline molybdenum disulfide (MoS₂), Modell. Simul. Mater. Sci. Eng., 2013, 21(4), 045003.
63. K. Jeong Won and H. Ho Jung, Comparison of C60 encapsulations into carbon and boron nitride nanotubes, J. Phys.: Condens. Matter, 2004, 16(23), 3901.
64. E. Bitzek, P. Koskinen and F. Gähler, et al., Structural Relaxation Made Simple, Phys. Rev. Lett., 2006, 97(17), 170201.
65. J. C. Zhang, D. C. Hou and R. M. Liu, et al., Effect of stacking order on the vibration properties of bilayer black phosphorus, Proc. R. Soc. A, 2022, 478(2263), 20220294.
66. D. Çakir, C. Sevik and F. M. Peeters, Significant effect of stacking on the electronic and optical properties of few-layer black phosphorus, Phys. Rev. B, 2015, 92(16), 165406.
67. Y. L. Mao, C. S. Xu and J. M. Yuan, et al., Effect of stacking order and in-plane strain on the electronic properties of bilayer GeSe, Phys. Chem. Chem. Phys., 2018, 20(10), 6929–6935.
68. Y. Liu, Y. Huang and X. F. Duan, van der Waals integration before and beyond two-dimensional materials, Nature, 2019, 567(7748), 323–333.

---