Tunable Negative Poisson's Ratio in Van der Waals Superlattice

Negative Poisson's ratio (NPR) materials are functional and mechanical metamaterials that shrink (expand) longitudinally after being compressed (stretched) laterally. By using first-principles calculations, we found that Poisson's ratio can be tuned from near zero to negative by different stacking modes in van der Waals (vdW) graphene/hexagonal boron nitride (G/h-BN) superlattice. We attribute the NPR effect to the interaction of pz orbitals between the interfacial layers. Furthermore, a parameter calculated by analyzing the electronic band structure, namely, distance-dependent hopping integral, is used to describe the intensity of this interaction. We believe that this mechanism is not only applicable to G/h-BN superlattice but can also explain and predict the NPR effect in other vdW layered superlattices. Therefore, the NPR phenomenon, which was relatively rare in 3D and 2D materials, can be realized in the vdW superlattices by different stacking orders. The combinations of tunable NPRs with the excellent electrical/optical properties of 2D vdW superlattices will pave a novel avenue to a wide range of multifunctional applications.


Introduction
Negative Poisson's ratio (NPR) material shrinks laterally when axially compressed or laterally expands when subjected to axial stretching. Compared to positive Poisson's ratio (PPR) materials, NPR material has higher indentation resistance [1], larger impact resistance, more superior sound absorption performance, and better crack propagation resistance [2]. In addition, when subjected to out-of-plane bending moments, the NPR material will exhibit a dome shape rather than the PPR material tending to saddle shape. These excellent properties indicate that the NPR materials have broad application prospects in the automotive, aerospace, marine, and other industrial fields [3].
Moreover, although many in-depth studies have examined the possible existence of NPR effects in 3D and 2D materials, few studies have reported the NPR effect in vdW superlattice. With the development of manufacturing technology, graphene-based superlattices showed enhanced stability in experiments [35]. Therefore, the designability and diversity of vdW superlattices provide a broad prospect for designing multifunctional NPR materials, such as NPR electrodes and molecular sieves [36]. In previous studies, the most NPR phenomena can be attributed to its unique wrinkle or re-entrant structures. In addition to these special geometric reasons, we want to study the fundamental physical mechanisms that form NPR effect.
In our previous study, we reported near-zero Poisson's ratio (ZPR) phenomena in G/h-BN and multilayer h-BN [37]. Interestingly, in this study, using first-principles calculations, we found that Poisson's ratios of G/h-BN superlattices are -0.109, -0.111, and 0.023 in different stacking modes. The dichotomy between NPR and PPR effects exhibited in the G/h-BN superlattice, which can be explained by a special electronic structure at the interfacial layer. Although G/h-BN is a kind of simple vdW heterostructure, it is convenient to make theoretical analysis and calculation clearly. It may open a beginning for the theoretical study of NPR effect in different stacking modes of vdW materials. In addition, we calculated the out-of-plane stiffness of the G/h-BN superlattice with different stacking modes. These modes with NPR also have outof-plane negative shear modulus (NSM), i.e., when shear strain is applied to NSM materials, as the shear strain increases, the corresponding shear stress tends to decrease.
Ultimately, the NPR phenomenon, which was relatively rare in 3D and 2D materials, can be realized in the vdW superlattices by different stacking orders as designed. Furthermore, studying on how to change the PPR material into the NPR material not only has important practical engineering application value but also theoretical value for in-depth study on other possible related interesting physical properties, such as negative pressure electricity, negative stiffness, and negative thermal expansion.

Results
The unit cell of the G/h-BN superlattice is composed of 1 × 1 graphene unit cell (2 C atoms) and 1 × 1 h-BN unit cell (1 B atom and 1 N atom) in the x-y plane. The lattice constant a of the pristine monolayer graphene and h-BN are 2.465 Å and 2.509 Å, respectively. Then, the a of G/h-BN superlattice calculated by first-principles calculations is 2.485 Å, and the lattice mismatch between graphene and h-BN is less than 1%.
The interlayer binding energies (E bind ) and equilibrium distances (d) of all stacking modes of G/h-BN superlattices obtained by density functional theory (DFT) approach are given in Table S1. Here, we investigated three highly symmetric stacking modes of G/h-BN superlattices: N atom sublattice on hexagonal C atom ring (stacking mode A), B atom sublattice on C atom ring center (stacking mode B), and N atom sublattice on C atom ring center (stacking mode C) [38]. E bind follows the order of E bind ðAÞ < E bind ðBÞ < E bind ðEÞ , while d follows the order of d A > d B > d C .
2.1. Stiffness. By analyzing the strain energy, the elastic constants of G/h-BN superlattices were derived from the linear fitting of the energy-strain relationship (Table S2). For hexagonal crystal, the in-plane mechanical properties of G/h-BN superlattice are isotropic (Y 11 = Y 22 , v 12 = v 21 , v 13 = v 23 ) [41]. Young's modulus Y αα is determined by elastic constants C αβ (see method section). Notably, the in-plane Y 11 of the 2D material is the product of the Y 11 of the corresponding 3D material and the effective thickness [27], and we took the d as the effective thickness for each layer of 2D material. Table 1 shows that the out-of-plane Y 33 of the stacking mode A, B, and C is 44.9, 45.6, and 49.0 GPa, respectively. The smaller the d of the stacking mode, the greater the corresponding Y 33 . In addition, we calculated the in-plane Y 11 of the G/h-BN superlattice, which is almost equal to the sum of Y 11 of the monolayer graphene and h-BN. Therefore, this result explains the reason for the stability enhancement of carbon-based superlattices observed in experiments [35]. However, the difference in Y 11 of the superlattice is mainly due to the different in d.

Poisson's Ratio.
We compared G/h-BN superlattices in different stacking modes under different uniaxial strains along the x direction (ε x ) ( Figure 1). For stacking modes A and B, the d is auxetic for ε x > 0, but the same phenomenon was not found in stacking mode C. Interestingly, stacking modes A and B not only have the NPR (v 13 ) effect but they also have negative shear modulus (NSM) (G 44 ) in the outof-plane direction. The shear force decreases with the increase of shear deformation, which is the NSM effect.
In order to study the anisotropy of Poisson's ratio of these materials, orientation-dependent Poisson's ratio was calculated ( Figure 2). We found that the stacking modes A, B, and C have the ZPR (v 13 ) effect at θ = 36:4°, 36.8°, and 18.7°, respectively. Therefore, stacking modes A and B exhibit a NPR effect in a larger crystal orientation angle range than stacking mode C.
To our knowledge, it is very difficult to measure Poisson' s ratio of several layers of two-dimensional (2D) material with the existing experimental method. Because for these ultra-thin films, when the in-plane strain is applied, the out-of-plane deformation is very small and difficult to observe. However, for multilayer 2D materials, X-ray diffraction can be easily used to measure Poisson' s ratio when the thickness is close to 10 nm [42]. The NPR effect is generated in the interfacial layer. Therefore, both multi and single-layer vdW materials can exhibit the same NPR. It is relatively easy to measure Poisson' s ratio for the multilayer vdW materials with a certain thickness.

Interlayer Binding Energy.
Assuming that the interaction between the two layers of the superlattice is additive, the binding potential can be expressed as the cumulative interaction of atoms between different layers [43]. The binding energy of two atoms combined by vdW forces can be expressed by the Lennard-Jones potential: Here, r represents the distance between the two atoms.
The ε and σ are fitting constants. The first term represents the vdW attraction, and the second term represents Pauli's For monolayer materials, the unit of Young's and shear modulus is Nm -1 . For G/h-BN superlattices, the unit of Young's and shear modulus is GPa.

Side view
Top view   repulsion [44]. Therefore, the interlayer potential of the vdW superlattice can be expressed as where ρ 1 and ρ 2 are the mass densities of two layers of vdW superlattice, respectively. The distance rðxÞ = ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi is obtained from the geometric relationship between the coordinate x and d. For G/h-BN superlattice, ρ 1 = ρ 2 . In Figure 2, the fitting curves of~d −4 below the horizontal coordinate axis represent vdW attraction, and the curves of~d −10 above the horizontal coordinate axis represent Pauli repulsion. Therefore, the low-order vdW term plays a major role in the large d, while high-order Pauli's repulsion term plays a major role in the small d.
According to the first-principles calculation, Eq. (2) can well describe the E bind of the vdW superlattice.  Figure 2(f)). The binding energy increases slightly; so, the G/h-BN superlattice exhibits the ZPR effect. According to the first-principles calculation, Poisson's ratio of the material has a relationship with σ. If the material can exhibit the NPR effect, the value of σ under tensile strain (σ′) is greater than the initial value of σ (Table S4). However, we calculated that the increase of ε is not a necessary condition for the NPR effect, i.e., σ plays a major role in the NPR effect.
3.2. Relationship between p z Orbitals and NPR. The Bloch wave function of the p z orbital electron in a periodic lattice under the tight binding (TB) approximation can be expressed as For G/h-BN superlattice, C atoms in graphene and N atoms in h-BN have p z orbitals. When the atom A is used as the origin of coordinates (Figures 3(a) and 3(b)), let the in-plane strain be the perturbation δ The constant C is a normalization constant, which should satisfy the normalization condition hψ p z ð k * Þ | ψ p z ð k * Þi = 1. The density-weighted length of the p z electrons in the outof-plane direction can be expressed as Here, l z ð k * Þ is the length of the p z electrons with momentum k * in the out-of-plane direction. Therefore, the length (L z ) of p z electrons with all momentum should be the integral of l z ð k * Þ in the first Brillouin zone (BZ). Finally, the charge density-weighted length of the p z orbital in the out-of-plane direction can be obtained: Þi, which is the length of the isolated p z orbital (Figure 3(c)). Therefore, according to Eq. (6), we got the analytical solution of the relationship between ε x and L z ðδ * Þ (Figure 3(d)). The calculation details can be found in the Supporting Information. Meanwhile, partially differentiate f ðδ * Þ to the in-plane perturbation δ * j is ∂f ðδ * Þ/∂ δ * > 0. Therefore, p z orbitals will extend out-of-plane under inplane tensile strain.
Since the length of the p z orbital has auxetic effect under in-plane strain, we quantitatively studied the charge density distribution along the out-of-plane direction of G/h-BN superlattice by using first-principles calculations. In the out-of-plane direction, the charge density at coordinate z can be expressed as Here, ρ E ðx, y, zÞ is the charge density at the coordinate (x, y, z) with the energy of E, and ε F is the Fermi level of the system. In order to quantify the change of the charge density in the out-of-plane direction under stress, we calculated the weighted length of the electron density in the out-of-plane direction according to the following formula: Notably, the charge of the graphene in the out-of-plane direction is mainly contributed by the p z orbital. The L z of each layer in the G/h-BN superlattice under in-plane strain ε x = 0 and ε x = 0:08 is shown in Table S4. In each layer of the G/h-BN superlattice, the L z is elongated. When an in-plane strain ε x = 0:08 was applied, the bond angle ∠NBN increased from 120°to 122.36°, resulting in the charge localization (Fig. S2). Quantitatively, we found that the value of L Z of the monolayer h-BN and graphene in G/h-BN superlattices increased by 1.8~1.9% and 2.3~2.4%, as the in-plane tensile strain increases by 8% (Table S5), explaining the NPR effect in stacking modes A and B along the out-of-plane direction. This is consistent with the analytical solution obtained by TB approximation (Figure 3(d)). For stacking mode C, the N atom sublattice is on the C atom ring center. In h-BN, the N atom has a fully filled p z orbital, while the B atom has an empty p z orbital. The p z orbitals of G/h -BN superlattices hardly overlap (Figure 4(c)); so, the change of the p z orbitals has little effect on the Pauli repulsion between the interfacial layers, resulting in no significant NPR effect. Figure 5 shows the DFT and TB-based band structure of G/h-BN superlattices in different stacking modes. To further understand the first-principles calculation results, we adopted the TB model to describe the electrons in G/h-BN superlattices with different stacking modes near the Fermi level. During the TB calculation, a unit cell contains two C 1 and C 2 carbon atoms at different positions and one N atom. Since the electronic states of the three bands around the Fermi level are completely contributed by the p z orbitals of C 1 , C 2 , and N atoms, only the p z orbitals of C 1 , C 2 , and N atoms are included in the TB model.

Relationship between Electronic Band Structures and NPR.
where subscripts 1, 2, and 3 represent C 1 , C 2 , and N atoms, respectively. Because the interlayer distance is longer than the C-C bond length, the nearest-neighbor interaction between C and N atoms and the next nearest-neighbor interaction between C and C atoms are considered (detailed Hamiltonian matrix elements can refer to the Supporting  Information). The distance-dependent hopping integral is determined by the formula

A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A B B B B B B B'
Here, d 0 represents the interfacial layer equilibrium distance, and V ijσ is the hopping integral between p z orbitals at d 0 . d ij is the distance between the ith and jth atoms, and q ij is the decay constant for the integral [45]. For the G/h-BN superlattices, the values of V ijσ and q ij can refer to the Supporting Information Table S6.
The distance-dependent hopping integral (f CNσ ) describes the intensity of the interaction between the p z orbitals of C and N atoms. Therefore, f CNσ is a power-exponential function of the interlayer spacing and is proportional to the NPR (Table 2). Furthermore, the higher the value of f CNσ , the greater the value of the corresponding NPR. Note that after Taylor expansion of the f CNσ ðd ij Þ, the quadratic term is the previous research results [46,47].
Consequently, the vdW superlattice can exhibit an NPR effect only if they have p z orbitals in the out-of-plane direction, and the p z orbitals overlap between the interfacial layers. Meanwhile, the NPR effect in all vdW materials can be explained by the same physical mechanism given in this  section. For example, for lattice-matched materials, a previous study showed that AA-stacked h-BN (a N atom on a N atom in another layer) can exhibit an NPR effect [26] ( Table 3). In addition, for lattice-mismatched vdW materials, the p z orbitals between the interfacial layers overlap; so, these materials should exhibit an NPR effect. For example, it has been observed that WS 2 /WSe 2 heterostructure expands abnormally under engineering tensile strain [48]. Therefore, according to this physical mechanism, the NPR phenomenon should exist in a large number of vdW materials, which was considered as a rare phenomenon in bulk and monolayer 2D materials.
Moreover, in the experiment, the isolated atomic layers can also be reassembled into the designed heterostructure layer by layer in a precisely selected order [49]. Therefore, for the same kind of investigated material, it can also be switched in different stacking modes through experimental methods. For example, the G/h-BN superlattices may also be tuned among stacking modes A, B, and C. Similarly, we can change the material without strong interlayer p z orbital interaction into a material with p z orbitals strongly overlapping between the interfacial layers, thus exhibiting an NPR effect.
In conclusion, we studied Poisson's ratios and the binding energies of G/h-BN superlattices in different stacking modes by using the first-principles method. We found that the stacking mode C has a ZPR effect at the interfacial layer, while the stacking modes A and B show NPR effects. The NPR effect is mainly due to the interaction of the p z orbitals between the interfacial layers. Furthermore, the distancedependent hopping integral (f CNσ ) calculated by analyzing that the electronic band structure can be used to describe the intensity of this interaction. The f CNσ is a powerexponential function of the interlayer spacing and is proportional to the NPR. Moreover, we calculated their Young's and shear modulus and found that the stacking modes A and B also have NSM effect in the out-of-plane direction. These materials with negative index coexistence will provide broad prospects for multifunctional and multipurpose materials. Finally, we expect that the theory can be verified by experiments and provide a solid foundation for the large-scale searching and predicting NPR materials in the future.

Methods
Based on density of functional theory, all first-principles calculations were implemented by the planewave projector augmented wave (PAW) method in Vienna ab initio simulation package (VASP) code [50]. The exchange correlation functional adopted the generalized gradient approximation (GGA) of the Perdew−Burke−Ernzerhof (PBE) functional [51]. In order to test the robustness of our results, the vdW-corrected functionals proposed by Grimme DFT + D2 [52], DFT + D3 [53], many-body dispersion (MBD) [54], and vdW-corrected functional optB88-vdW [55] methods were used in first-principles calculations. In this paper, the calculation results of functional optB88-vdW are given because of its good agreement with the experimental results, and the results obtained by different vdW-corrected methods are only slightly different in numerical value (the detailed results are in the Supporting Information).
The G/h-BN superlattice was calculated by using 28 × 28 × 10 Monkhorst-Pack K-point mesh. The energy cut-off value is 500 eV, and the structures were completely relaxed until their atomic Hellmann-Feynman forces were less than 0.005 eV/Å. The convergence criterion of energy in the selfconsistency process is 10 -6 eV. We also calculated electronic band structures for G/h-BN superlattices by using the HSE06 hybrid functional [56].
To quantitatively characterize the mechanical properties of the interface, the interlayer binding energy (E bind ) between the monolayer graphene and h-BN is as follows: where E G/h−BN , E G , and E h−BN are the energies of the G/h-BN superlattice, graphene, and h-BN, respectively. S represents the in-plane area of the superlattice. The elastic constant is defined by expanding the internal energy E into Taylor series in elastic strain at constant entropy. The expansion coefficient in the Taylor series is the elastic constant [57]: where ρ 0 and η ij are the initial mass density and the Lagrangian strains of the material [58]. In this work, we use contracted notations (11 → 1, 22 → 2, 33 → 3, 13 → 4, 23 → 5, 12 → 6, C ijkl → C αβ ) as tensor indices. In addition, we Young's modulus for the material is computed by Poisson's ratio is defined as where ε i is the strain in the direction of uniaxial loading (in the i-direction), and ε j is the resulting strain in the transverse direction (the j-direction). In our calculations, we applied different uniaxial strains to the lattice. This strained structure was then completely relaxed to evaluate the magnitude of the strain in the out-of-plane direction. The detailed calculation process of the relationship between θ and ν 13 is provided in the Supporting Information.

Data Availability
All data needed to evaluate the conclusions in the paper are presented in the paper and supplementary materials. And additional data are available from the corresponding authors upon reasonable request.

Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.