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 p z 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.

1. 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].

The NPR phenomena have been found in many natural and artificial materials, such as cubic metals [4], porous polymers [5], honeycombs [6], perovskite [7, 8], silicon oxides [7, 8], ceramic aerogels [9], reentrant crystal structures [10-17], and liquid crystal elastomer [18]. Recently, the NPR effect has also been found in some two-dimensional (2D) materials [19], such as black phosphorus [20, 21], δ-phosphorene [22], borophene [23, 24], graphene [25], h-BN [26], 1 T-type transition metal dichalcogenides [27], group-IV monochalcogenides [28], Be5C2 [29], silicon dioxide [30, 31], FeB6 [32], B4N [33], and Ag2S [34]. In addition, there are hundreds of thousands of materials in the inorganic crystal structure database (ICSD); however, the number of NPR materials reported in the study is a few hundred, which is relatively small. Among them, artificial materials and structures often have a very high NPR, while crystal materials have a relatively small NPR. Especially in 2D materials, NPR is smaller. In 2014, NPR was first found in 2D materials, namely, black phosphorus, with a value of -0.027 [20]. After that, the NPR of the 2D materials found in the later study mostly remains near -0.1. Based on our theoretical calculations, unlike pure 2D materials, the NPR in the vdW heterostructure can only be maintained around -0.1, which is due to the expansion amplitude of p orbital under in-plane strain.

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 out-of-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.

2. 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 (Ebind) 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]. Ebind follows the order of Ebind (A) < Ebind (B) < Ebind (E), while d follows the order of dA > dB > dC.

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 (Y11 = Y22, v12 = v21, v13 = v23) [41]. Young's modulus Y is determined by elastic constants C (see method section). Notably, the in-plane Y11 of the 2D material is the product of the Y11 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 Y33 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 Y33. In addition, we calculated the in-plane Y11 of the G/h-BN superlattice, which is almost equal to the sum of Y11 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 Y11 of the superlattice is mainly due to the different in d.

Table 1

Poisson's ratio, Young's, and shear modulus of graphene, h-BN and G/h-BN.

		Method	In-plane			Out-of-plane
			v 12	Y 11	G 66	v 13	v 31	Y 33	G 44
G/h-BN	Stacking mode A	DFT	0.199	886.8	369.7	-0.109	-0.005	44.9	-1.3
	Stacking mode B	DFT	0.198	896.4	374.2	-0.111	-0.006	45.6	-2.8
	Stacking mode C	DFT	0.199	932.8	389.1	-0.023	-0.001	49.0	8.7
Monolayer	Graphene	Expt. [39]	0.165	340 ± 50	145.9 ± 30
		DFT	0.159	340	146.7
	h-BN	Expt. [40]	0.19	273	114.7
		DFT	0.199	238	99.3

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.

2.2. Poisson's Ratio

We compared G/h-BN superlattices in different stacking modes under different uniaxial strains along the x direction (ε) (Figure 1). For stacking modes A and B, the d is auxetic for ε > 0, but the same phenomenon was not found in stacking mode C. Interestingly, stacking modes A and B not only have the NPR (v13) effect but they also have negative shear modulus (NSM) (G44) in the out-of-plane direction. The shear force decreases with the increase of shear deformation, which is the NSM effect.

Figure 1

Top and side views of G/h-BN superlattices in different stacking modes under in-plane strain (a, d, g) ε = −0.08, (b, e, h) ε = 0, and (c, f, i) ε = 0.08. Here, d represents interfacial layer equilibrium distance.

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 (v13) 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.

Figure 2

Poisson's ratio v(θ) (v13) as a function of G/h-BN superlattices in (a) stacking mode A, (b) stacking mode B, and (c) stacking mode C for i fixed in the x direction and j varying in the y-z plane. Interlayer binding energy (Ebind) of G/h-BN superlattices in (d) stacking mode A, (e) stacking mode B, and (f) stacking mode C with d. The fitting energy curves below the horizontal coordinate axis represent vdW attraction, and the fitting energy curves above the horizontal coordinate axis represent Pauli repulsion. These dashed lines represent the fitting energy curves for ε = 0.08, and the solid lines represent the fitting energy curves at the equilibrium position. The insets show enlarged energy curves.

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.

3. Discussion

3.1. 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 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 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 Ebind of the vdW superlattice.

Figures 2(d) and 2(e) show that when ε = 0.08, the Pauli repulsion energy increases significantly, while the vdW attraction energy has a negligible change. Therefore, with the decrease of Ebind, d of the G/h-BN superlattice expands, while the lowest point of the energy curve moves forward along the positive direction of the x axis, resulting in a negative v13 of -0.109 and -0.111. For the stacking mode C, vdW and Pauli repulsion show negligible changes under ε = 0.08 (see 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 Orbitals and NPR

The Bloch wave function of the p 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 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 wave function of p orbitals of the atom A is . Meanwhile, the wave function of p orbitals of the atom B is . When the A and B atoms in the lattice are bonded, the wave function of the bonded p orbital can be expressed by the linear combination of atomic orbitals as The constant C is a normalization constant, which should satisfy the normalization condition . The density-weighted length of the p electrons in the out-of-plane direction can be expressed as Here, is the length of the p electrons with momentum in the out-of-plane direction. Therefore, the length (L) of p electrons with all momentum should be the integral of in the first Brillouin zone (BZ). Finally, the charge density-weighted length of the p orbital in the out-of-plane direction can be obtained: where and , which is the length of the isolated p orbital (Figure 3(c)). Therefore, according to Eq. (6), we got the analytical solution of the relationship between ε and (Figure 3(d)). The calculation details can be found in the Supporting Information. Meanwhile, partially differentiate to the in-plane perturbation is . Therefore, p orbitals will extend out-of-plane under in-plane tensile strain.

Figure 3

Unit cell of (a) graphene and (b) h-BN. (c) Wave function superposition of A and B atoms. (d) The length of p orbital in each layer of G/h-BN superlattice under different strains. Here, a represents the distance between A and B atoms. Under the in-plane strain ε = δ/a, the length of p orbital is L + δL.

Since the length of the p 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, ρ(x, y, z) is the charge density at the coordinate (x, y, z) with the energy of E, and ε 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 orbital. The L of each layer in the G/h-BN superlattice under in-plane strain ε = 0 and ε = 0.08 is shown in Table S4. In each layer of the G/h-BN superlattice, the L is elongated. When an in-plane strain ε = 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 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 orbital, while the B atom has an empty p orbital. The p orbitals of G/h-BN superlattices hardly overlap (Figure 4(c)); so, the change of the p orbitals has little effect on the Pauli repulsion between the interfacial layers, resulting in no significant NPR effect.

Figure 4

The charge density with the isosurface of 0.103 Å−3 and the electron localization function (ELF) with the Miller indices of (1 1 0) of h-BN superlattices in (a) stacking mode A, (b) stacking mode B, and (c) stacking mode C. The red dotted arrows and circles indicate the overlap of the p orbits, while the black dotted arrows and circles show that the overlap of the p orbitals does not actually exist.

3.3. Relationship between Electronic Band Structures and NPR

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 C1 and C2 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 orbitals of C1, C2, and N atoms, only the p orbitals of C1, C2, and N atoms are included in the TB model. The Hamiltonian matrix can be written as where subscripts 1, 2, and 3 represent C1, C2, 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 Here, d0 represents the interfacial layer equilibrium distance, and V is the hopping integral between p orbitals at d0. d is the distance between the ith and jth atoms, and q is the decay constant for the integral [45]. For the G/h-BN superlattices, the values of V and q can refer to the Supporting Information Table S6.

Figure 5

Electronic band structures of G/h-BN superlattices in different stacking modes from DFT (gray) and TB (black) calculations. Red, green, and blue denote the contributions of the p orbitals of C1, C2, and N atoms, respectively.

The distance-dependent hopping integral (f) describes the intensity of the interaction between the p orbitals of C and N atoms. Therefore, f is a power-exponential function of the interlayer spacing and is proportional to the NPR (Table 2). Furthermore, the higher the value of f, the greater the value of the corresponding NPR. Note that after Taylor expansion of the f(d), the quadratic term is the previous research results [46, 47].

Table 2

The values of d0, f, and Poisson's ratio of G/h-BN superlattices in different stacking modes.

Stacking mode	d 0 (Å)	f CNσ (eV)	Poisson's ratio
A	3.448	-0.31	-0.109
B	3.416	-0.31	-0.111
C	3.274	-0.22	-0.023

Consequently, the vdW superlattice can exhibit an NPR effect only if they have p orbitals in the out-of-plane direction, and the p 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 orbitals between the interfacial layers overlap; so, these materials should exhibit an NPR effect. For example, it has been observed that WS2/WSe2 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.

Table 3

Poisson's ratios of vdW materials with NPR.

	v 13
Bilayer graphene [26]	-0.09
AA-stacked h-BN [26]	-0.12
G/MoS2 heterostructure [37]	-0.09
G/h-BN
Stacking mode A	-0.109
Stacking mode B	-0.111

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 orbital interaction into a material with p 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 orbitals between the interfacial layers. Furthermore, the distance-dependent hopping integral (f) calculated by analyzing that the electronic band structure can be used to describe the intensity of this interaction. The f is a power-exponential 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.

4. 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 self-consistency 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 (Ebind) between the monolayer graphene and h-BN is as follows: where EG/, EG, and E 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 η 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 → C) as tensor indices. In addition, we define Lagrangian strains η1 → ε, η2 → ε and η3 → ε. The compliance coefficients S are defined as Young's modulus for the material is computed by Poisson's ratio is defined as where ε is the strain in the direction of uniaxial loading (in the i-direction), and ε 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.