Negative Poisson’s ratio in two-dimensional honeycomb structures

Negative Poisson’s ratio (NPR) in auxetic materials is of great interest due to the typically enhanced mechanical properties, which enables plenty of novel applications. In this paper, by employing first-principles calculations, we report the emergence of NPR in a class of two-dimensional honeycomb structures (graphene, silicene, h-BN, h-GaN, h-SiC, and h-BAs), which are distinct from all other known auxetic materials. They share the same mechanism for the emerged NPR despite the different chemical composition, which lies in the increased bond angle (θ). However, the increase of θ is quite intriguing and anomalous, which cannot be explained in the traditional point of view of the geometry structure and mechanical response, for example, in the framework of classical molecular dynamics simulations based on empirical potential. We attribute the counterintuitive increase of θ and the emerged NPR fundamentally to the strain-modulated electronic orbital coupling and hybridization. It is proposed that the NPR phenomenon can also emerge in other nanostructures or nanomaterials with similar honeycomb structure. The physical origin as revealed in our study deepens the understanding on the NPR and would shed light on future design of modern nanoscale electromechanical devices with special functions based on auxetic nanomaterials and nanostructures.


INTRODUCTION
For common materials, the lattice along one direction expands or shrinks as the other orthogonal directions are compressed or stretched, respectively. The Poisson's ratio (ν) is the physical property used to quantify the phenomena, which is a positive value for most cases. However, for some special materials (the socalled auxetic materials), the ν can be negative 1 . The negative Poisson's ratio (NPR) is of great interest because the type of materials with NPR typically possess enhanced toughness, shear resistance, and efficient sound/vibration absorption, which enables plenty of novel applications, such as aerospace and defense 2 . Historically, researches of the NPR are mostly conducted on the bulk auxetic structures, which are further extended to nanomaterials and more interesting phenomena are found. For instance, NPR is recently found in carbon nanotubes and metal nanoplates 3,4 . In addition, since the discovery of graphene, twodimensional (2D) materials have been intensively studied for many years with promising applications in various fields 5 . The NPR has been also realized in 2D materials with specific engineering, such as by cutting into nanoribbons 6 , introducing vacancy defects 7 , creating periodic porous 8 , rippling curvature at very high temperatures 9 , etc. Moreover, the NPR has been also reported recently in 2D materials without any external modifications to the structure, shape, or composition, such as the in-plane NPR (for instance, borophene 10 , penta-graphene 11 , α-silica 12 , Be 5 C 2 13 ), the out-of-plane NPR (phosphorene 14,15 , GeS 16 , SnSe 17 , arsenic 18,19 , TiN 20 ), and both in-plane and out-of-plane NPR possessed in bidirectional 2D auxetic material Ag 2 S 21 .
Currently, the understanding of the NPR phenomena is dominated by the geometry analysis in literature 1,4,6,7,14,[22][23][24] . The auxetic effect is generally thought to be independent of chemical composition and electronic structure, which originates from the special reentrant structures or the rigid building blocks linked by flexible hinges. However, Yu et al. 25 reported recently that in-plane NPR emerges in the 1T-type 2D transition-metal dichalcogenides (TMDCs), where the crystal structure possesses no reentrant or hinge-like building block. The electronic effect rather than the mechanical factor is found responsible for the NPR 25,26 . Therefore, an interesting and thought-provoking question arises that are there other exceptional nanomaterials without reentrant or hinge-like building block possessing the NPR beyond the TMDCs?
In this paper, based on first-principles calculations, we systematically investigate the mechanical response to uniaxial strain of a class of 2D honeycomb structures (graphene, silicene, h-BN, h-GaN, h-SiC, and h-BAs). By applying uniaxial strain along the zigzag and armchair directions, respectively, it is found the NPR emerges in all the six 2D materials, but only with the strain applied along the armchair direction. The honeycomb structure with emerged NPR is distinct from all other known auxetic materials. By analyzing the evolution of key geometry parameters (bond length and angle) with strain increasing, it is found that the emerged NPR is due to the anomalous increase of bond angle, which, however, cannot be explained in the traditional point of view of the geometry structure and mechanical response. Deep understanding on the counterintuitive increase of bond angle is achieved from a fundamental view of electronic structure based on the analysis of orbital pDOS beyond the mechanical analysis. Our study on the NPR in the class of 2D honeycomb materials deepens the understanding on the origins of NPR, which would shed light on the advanced design of modern nanoscale electromechanical devices. 1 Figure 1 shows the response of driven strain, stress, and energy change per atom to the applied strains for the six 2D materials of graphene, silicene, h-BN, h-GaN, h-SiC, and h-BAs. When the uniaxial strains are applied along the zigzag and armchair directions, respectively, the mechanical response is anisotropic for all the six 2D materials. The stress is scaled by substituting the thickness including vacuum space with the effective layer thickness (graphene: 3.5 Å; silicene: 4.65 Å; h-BN: 3.1 Å; h-GaN: 3.74 Å; h-SiC: 4.2 Å; h-BAs:3.7 Å) based on the atomic van der Waals diameters [27][28][29][30][31] .

Mechanical properties
Generally, the lattice constant along one direction decreases when the tensile strain is applied along other orthogonal directions. However, in graphene, the lattice constants along both the zigzag and the armchair directions simultaneously increase with large uniaxial tensile strains applied along either the zigzag (Fig. 1a) or the armchair (Fig. 1d) directions. While for silicene, h-BN, h-GaN, h-SiC, and h-BAs, the simultaneous increase of the lattice constants is only found when large uniaxial tensile strains are applied along the armchair direction, as shown in Fig. 1d. The simultaneous increase of the lattice constants with uniaxial strain applied is anomalous for its derived NPR. Note that the simultaneous increase of the lattice constants for silicene is very small compared to graphene, h-BN, h-GaN, h-SiC, and h-BAs, which might be due to the effect of its buckling structure.
With much larger strains applied, the structures fail, which is revealed by the mutation of lattice constants (Fig. 1a, d). The corresponding abrupt decreases of stress (Fig. 1b, e) and energy change per atom (Fig. 1c, f) also reveal the failure of the structure at very large strains, and the failure points are consistent with each other, respectively (see Supplementary Note 1 for more information). The situation for h-GaN, h-SiC, and h-BAs is special compared to the other three 2D materials that their structures never fail when the uniaxial strain is applied along the armchair direction, which might be due to the strongly polarized bonds 29,31,32 .
Negative Poisson's ratio The Poisson's ratio is defined as 22 where ϵ 1 is the driven strain due to the applied strain (ϵ 2 ). Since remarkable NPR only arise with uniaxial strain applied along the armchair direction as indicated in Fig. 1a, d, we extract the Poisson's ratio in Fig. 2 for the six 2D materials with strains applied along the armchair direction based on the results as shown in Fig.  1d. Certainly, we only consider the cases with strain smaller than the failure points. It shows that NPR emerges when the applied uniaxial strain is larger than some threshold values, where the lateral lattice constants begin to increase. The threshold values (~18%) are close to each other for graphene, h-BN, h-GaN, h-SiC, and h-BAs, despite their different NPR values and varying trends. In contrast, the threshold value is much larger for silicene (33%) and the NPR value is much smaller, which is consistent with the variation of lattice constant as analyzed above and could be attributed to the effect of its distinct buckled structure. Besides, the magnitude of NPR in h-GaN is unexpectedly large compared to graphene, silicene, h-BN, phosphorene 14 , and other lots of materials 1 .
Previous study from classical MD simulations reported that the NPR in graphene emerges when the strain along the armchair direction exceeds 6% 22 . The discrepancy of the thresholds might be due to the different computational methods employed (firstprinciples vs. classical MD with empirical potential).

Mechanistic explanation
Despite the different components, NPR emerges in all the six 2D materials with tensile strain applied along the armchair direction, which possess similar honeycomb structure. The honeycomb structure possesses no reentrant or hinge-like building block, which is distinct from all other known auxetic materials 1 . To understand how the NPR emerges, we study the evolution of the key geometry parameters to gain some insight. As shown in the inset of Fig. 3a, the lattice constant along the zigzag direction (l zigzag ) is governed by the bond length of b 1 and the bond angle of θ. Note that we only consider the projection in the xy-plane of b 1 and θ for silicene with buckling structure. With tensile strain applied along the armchair direction, the variation of bond length b 1 and bond angle θ shows similar behavior for the six 2D materials, as shown in Fig. 3. The bond length b 1 first increases and then decreases (Fig. 3a). While for the bond angle θ, it first Fig. 1 The response to strain. The anisotropic response of the a, d driven strain, b, e stress, and c, f energy change to the uniaxial strains applied along the a−c zigzag and d−f armchair directions, respectively. The simultaneous increase of driven strain and applied uniaxial strain is found for all the six 2D materials with honeycomb structure when the strain is applied along the armchair direction. Fig. 2 The negative Poisson's ratio. The calculated Poisson's ratio for the six 2D materials with honeycomb structure when the uniaxial tensile strain is applied along the armchair direction.
decreases and then increases, which is opposite to the variation of b 1 . The opposite variation of θ and b 1 can be intuitively understood in terms of their relation in geometry structure. When θ increases, the stretching force projected along the bond b 1 decreases, leading to the decreased bond length. As shown in the inset of Fig. 3a, the l zigzag is two times the projections of b 1 along the zigzag direction The variation of l zigzag is positively correlated to the variation of both b 1 and θ. However, the variation of l zigzag (Fig. 1d) is consistent with the variation of θ (Fig. 3b), while opposite to the variation of b 1 (Fig. 3a). Based on Eqs. (1) and (2), the relationship of Poisson's ratio with the variation of b 1 and θ is Thus, considering the large pristine magnitude of b 1 for the six 2D materials and the larger variation of θ than b 1 as shown in Fig.  3, the variation of l zigzag is dominated by the variation of θ, which would be slightly affected by the variation of b 1 . Consequently, the first decrease and then increase of the l zigzag lead to the emerged NPR when the strain along the armchair direction is larger than some threshold values.
It was found in a previous study based on classical MD simulations that 22 the bond stretching other than the angle bending is responsible for the NPR in graphene with strain applied along the armchair direction. The mechanism is in contrast to the conclusions as analyzed above in this study that the NPR in graphene is emerged due to the bond angle (θ) bending. The difference may lie in the difference between the classical MD simulations and the first-principles calculations. For instance, in classical MD simulations, the accuracy largely depends on the employed empirical potential that is used to describe interatomic interactions. While in first-principles calculations, the interatomic interactions are governed by the atomic pseudo-potential and the electronic structures. The empirical Brenner potential used in previous classical MD simulations does not involve any electronic properties. Thus, the widely used empirical Brenner potential might be not able to precisely describe the effect of strainmodulated electronic structures in graphene, and further the evolution of the geometry parameters (bond length b 1 and bond angle θ). Such limitation leads to the misinterpretation of the underlying mechanism for the NPR in graphene by the classical MD simulations, despite its success in capturing the NPR phenomenon.
In short, some insight have been achieved into the reason why the NPR emerges based on the analysis of the key geometry parameters evolutions (b 1 and θ) as the applied strains increase. With tensile strain applied along the armchair direction, the six 2D materials share the same mechanism for the emerged NPR despite the different components, which lies in the increased bond angle θ. However, the increasing θ with tensile strain applied along the armchair direction is counterintuitive as shown in Fig. 4h. The anomalous phenomena cannot be explained only in the traditional view of mechanical response of the geometry structure. Therefore, the insight from a fundamental view into the counterintuitive increasing θ (Fig. 4h) is necessary to gain deep understanding on the emerged NPR in the class of 2D materials, where the effect of the electronic structures must be involved.
Insight from electronic structure To reveal the underlying mechanism, we further study the evolution of orbital projected density of states (pDOS) of the six 2D materials with strains applied along the armchair direction. We plot the pDOS for comparison of the pristine state, the state where NPR is going to arise and the state where the structure is going to fail (Graphene: 0%, 14%, 27%; Silicene: 0%, 30%, 38%; h-BN: 0%, 13%, 25%; h-GaN: 0%, 17%, 33%; h-SiC: 0%, 16%, 37%; h-BAs: 0%, 26%, 35%). Generally for the six 2D materials, the contribution of p x orbital (along the zigzag direction) change nonmonotonically with the increasing strain applied along the armchair direction, which first increases slightly and then decreases largely, as shown in Fig. 4 for the three representative cases. The increased and decreased contribution are evidently shown by the approaching to and deviating away the valance band maximum (VBM), respectively. Consequently, the attractive interactions along the zigzag (x) direction are firstly enhanced and then weakened with the strain-modulated p x orbital coupling, as revealed by the orbitals' evaluation close to the VBM (Fig. 4a-m). For the 2D honeycomb structures studied in this work, all their p x -DOS near VBM present similar trend (first increase and then decrease), which reveals the changes of strain-modulated p x orbital coupling. Considering the geometry structure of the honeycomb structure as shown in Fig. 3a, the bond angle θ first decreases and then increases, leading to the emerged NPR.
It is worth to note the special cases of silicene and h-GaN. For silicene, due to the distinct buckling structure, there exists coupling between p x and p z orbitals, which is different from graphene and h-BN with planar structure. Consequently, the contribution of p x orbital is mediated by p z orbital, which does not decrease largely with strain applied. Thus, the interaction is not remarkably weakened, the increase of l zigzag in silicene is very small, and the NPR effect is very weak. As for h-GaN, previous studies have indicated the significant effect of Ga-d orbital 29,31 . The d orbital splits into two groups of d xy;yz;zx (t 2g ) and d x 2 Ày 2 ;z 2 e g À Á . Here we focus on the e g orbitals, which couple strongly with the p x orbitals. With the applied strain increasing, the coupling strength of e g -p x -orbitals becomes stronger due to the increased contribution of e g , which approaches the VBM as shown in Fig. 4f. The e gp x -orbitals coupling leads to the strongly polarized Ga-N bonding and protects the Ga-N bond from breaking. Thus, the structure of h-GaN never fails with the applied strain and the magnitude of NPR in h-GaN is unexpectedly large. Similar phenomena are also observed for h-SiC and h-BAs, which possess similar polarized bonding as h-GaN.
Based on the above discussion, the emergence of NPR in the six 2D materials together with the counterintuitive increase of θ are thoroughly understood from the point of view of electronic structure based on the analysis of orbital pDOS. Furthermore, the Fig. 3 The mechanical mechanism for the NPR. a, b The evolution with applied strain of the key geometry parameters (b 1 and θ) in the six 2D materials with honeycomb structure when the uniaxial tensile strain is applied along the armchair direction. The b 1 and θ are labeled on the unit cell in the inset of (a) for clarification.
NPR phenomenon can be anticipated to also emerge in other nanostructures or nanomaterials with similar honeycomb structure.

DISCUSSION
In summary, based on first-principles calculations, we performed systematic study on the anisotropic mechanical response to uniaxial strain of a class of 2D honeycomb structures (graphene, silicene, h-BN, h-GaN, h-SiC, and h-BAs). Intrinsic NPR is found in all the six 2D materials when the strain is applied along the armchair direction. The honeycomb structure is distinct from all other known auxetic materials. By analyzing the evolution of key geometry parameters (bond length (b 1 ) and bond angle (θ)) with strain increasing, it is found that the six 2D materials share the same mechanism for the emerged NPR despite the different components, which lies in the increased bond angle θ. However, the increase of bond angle θ at large strains is quite intriguing and anomalous, which cannot be explained in the traditional point of view of the geometry structure and mechanical response, such as in the framework of classical MD simulations based on empirical potential. Deep understanding on the counterintuitive increase of θ is achieved from a fundamental view of electronic structure based on the analysis of orbital pDOS. We find that the emerged NPR with the increased θ is fundamentally due to the strain-modulated orbital coupling and hybridization, which lead to the weakened attractive interactions. From the above results and conclusions, we propose that the NPR phenomenon can also emerge in other nanostructures or nanomaterials with similar honeycomb structure. Our study not only makes a comprehensive investigation of the NPR in the six 2D materials with honeycomb structure, but also reveals the physical origins, which deepens the understanding on the NPR and would shed light on future design of modern nanoscale electromechanical devices with special functions based on auxetic nanomaterials and nanostructures. Fig. 4 The electronic origin of the NPR. a-g, i-m The evolution of p x and e g orbital projected density of states (pDOS) with representative strains applied along the armchair direction (pristine → where NPR is going to arise → where structure is going to fail) for the six 2D materials with honeycomb structure. Inset of (b): The lateral view of the buckling structure of silicene. h A sketch for the anomalous increase of bond angle with stretch strain applied along the armchair direction.

METHODS
All the first-principles calculations are performed based on the density functional theory (DFT) using the projector augmented wave method 33 as implemented in the Vienna ab initio simulation package (VASP) 34 . The Perdew−Burke−Ernzerhof 35 of generalized gradient approximation is chosen as the exchange-correlation functional. The kinetic energy cutoff of wave functions is set as 2.5 times the maximal energy cutoff in the pseudopotentials for each material. The orthorhombic supercell is constructed containing two primitive cells (four atoms). The Monkhorst-Pack 36 k-mesh of 19 × 11 × 1 is used to sample the Brillouin Zone. The energy convergence threshold for the self-consistent field calculations is set as 10 −6 eV. The vacuum spacing of 20 Å is employed along the out-of-plane direction. The uniaxial strains are applied along the typical zigzag and armchair directions, respectively. The strain is defined as (l − l 0 )/l 0 , where l is the lattice constant under strains and l 0 is the original lattice constant with no strain. All geometries are fully optimized for all cases until the Hellmann−Feynman forces on all the atoms are smaller than 10 −5 eV/Å. By testing the effect of the size of the supercell in calculations, it is verified that the mechanical properties studied here are accordant. The structural stability of the strained situations where the NPR already emerges is verified by the 1000 fs DFT level MD simulations with the NVE ensemble (see Supplementary Materials for more information).