Auxetic two-dimensional transition metal selenides and halides

Auxetic two-dimensional (2D) materials provide a promising platform for biomedicine, sensors, and many other applications at the nanoscale. In this work, utilizing a hypothesis-based data-driven approache, we identify multiple materials with remarkable in-plane auxetic behavior in a family of buckled monolayer 2D materials. These materials are transition metal selenides and transition metal halides with the stoichiometry MX (M = V, Cr, Mn, Fe, Co, Cu, Zn, Ag, and X = Se, Cl, Br, I). First-principles calculations reveal that the desirable auxetic behavior of these 2D compounds originates from the interplay between the buckled 2D structure and the weak metal–metal interaction determined by their electronic structures. We observe that the Poisson’s ratio is sensitive to magnetic order and the amount of uniaxial stress applied. A transition from positive Poisson’s ratio (PPR) to negative Poisson’s ratio (NPR) for a subgroup of MX compounds under large uniaxial stress is predicted. The work provides a guideline for the future design of 2D auxetic materials at the nanoscale.


INTRODUCTION
As a fundamental mechanical property, Poisson's ratio (PR) describes the negative ratio of transverse strain to longitudinal strain. Most materials host positive Poisson's ratio (PPR), which means that they undergo a transverse contraction (expansion) when stretched (compressed). In contrast, materials with negative Poisson's ratio (NPR), i.e., the so called auxetic materials 1 , expand laterally when stretched and contract laterally when compressed. This unusual mechanical behavior has attracted significant attentions owing to its intriguing properties for both fundamental research and potential applications in biomedicine 2 , fasteners 3 , sensors 4 , national security and defense 5 , and other fields [6][7][8] .
Although the deformation mechanism of 2D materials with both positive and negative PR has been demonstrated in the literature, the intrinsic driving force, on the electronic structure level, for the auxetic behavior remains unknown in most cases. In our recent study 19 , it is proposed that the auxetic effect in group-VI 1T-type monolayer transition metal dichalcogenides (TMDs) is related to the gradually filling d-orbitals, which in turn determine the p-d bonding interactions. This observation indicates that the auxetic effect in 2D materials originates from the interplay of geometric and electronic structure effects. Unfortunately, although some metastable 1T-type structures have also been synthesized 21,22 , most TMD compounds that are predicted to host the NPR in the 1 T phase actually favor the 1H phase thermodynamically, which hinders experimental verification and demonstration.
In this article, based on the geometric evolution mechanism of 2D structures with a buckled square lattice, we propose a feasible approach to realize the NPR in a large class of 2D materials. Out of over 800 layered materials identified from the ICSD database 23 , motivated by a symmetry-based hypothesis, we identify a family of buckled square monolayer structures in the p4mm plane group that can host auxetic behaviors. Sharing the same stoichiometry of MX, these 2D materials are p4mm-type transition metal selenides and halides. Among them, monolayer compounds with d 10 cations tend to host in-plane NPR, whereas others usually exhibit PPR with small strain (<2%). This trend can be understood by analyzing the interactions between neighboring transition metal atoms up to the next-nearest neighboring (NNN) level. Experimental realization is called for this class of auxetic 2D compounds, as the bulk materials of several MX compounds considered in this work have been successfully synthesized [24][25][26][27] .

RESULTS AND DISCUSSION
Design rule of 2D NPR materials in a square lattice Our hypothesis is that the NPR effect in 2D materials can be achieved through a geometric evolution under strain in one of the simplest 2D geometric structures: the square lattice. Taking the planar p4mm-type structure as an example, as shown in Fig. 1a, a tensile strain along x-direction will give rise to an increasing distance between the two nearest neighboring (NN) atoms, thus inducing a downward (upward) resultant force for M 3 (M 4 ) atom. A transverse contraction along y-direction thus occurs, leading to a PPR. Obviously, a compressive strain along x-direction results in the upward (downward) movement of M 3 (M 4 ) atom, leading to a PPR as well. 1 Next, we illustrate that the NPR can originate from a more complicated structure by introducing an extra atom above/below the xy-plane, as shown in Fig. 1b, which has strong bonding interactions with its four neighboring atoms. Under a tensile strain along x-direction, the X atom moves downward owing to the attractive interaction with M 1 and M 2 atoms. If the NNN interactions between M atoms are negligible in comparison with the interaction between M and X atoms, the repulsive forces on M 3 and M 4 atoms from X atom dominates and an auxetic effect (NPR) thus occurs. However, if the NNN M-M interactions are not negligible, the M 3 and M 4 atoms would be subject to repulsive forces from the X atom and attractive forces from M 1 and M 2 atoms simultaneously. The resultant force induces a clockwise (counter-clockwise) rotation for the bond between M 3 (M 4 ) and X atom in the yz-plane that moves M 3 and M 4 atoms downward relative to M 1 and M 2 atoms. In this case, the tensile strain along xdirection induces a transverse contraction of the structure in xyplane (PPR). We expect that the electronic structure of host materials (such as the electronic states of metal atoms) plays a major role in the NNN interaction among M atoms, thus affecting the PRs of this set of 2D materials.
Poisson's ratios of 2D p4mm-type MX materials Based on the above hypothesis that the auxetic effect can be realized through the interplay of atomic geometry and electronic structure in a 2D buckled square lattice, we identify 2D compounds with square lattice and buckled geometric structures through a plane group symmetry-based data mining process from over 800 layered materials identified from the ICSD database 23 . The symmetry-based screening process is shown in Supplementary Figure 1. A family of compounds, sharing the MX stoichiometry (M: V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Ag; X: Se, Cl, Br, I) turn out to be the most promising candidates within the p4mm plane group. As shown in Fig. 2a, this family of 2D compounds exhibit buckled structures with a layer of metal cations (M: V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Ag) sandwiched between two layers of anions (X: Se, Cl, Br, I). Each cation M and its four neighboring anions X form a quasi-regular tetrahedron. Taking the X anion as a central point, a fourfold rotational symmetry is present in the xyplane and two mirrors are along x-and y-directions, which effectively place this class of 2D compounds in the p4mm plane group. The optimized lattice constants, band gaps, and elastic stiffness of the whole family obtained by HSE06 are illustrated as Supplementary Table 1.
The calculated PRs under small strain are shown in Fig. 2b. Note that the PRs along x-and the y-directions are equivalent, owing to the fourfold rotational symmetry. A striking observation is that the MX compounds with d 10 transition metal cations generally exhibit NPR, where AgCl has the most NPR at −0.18. The only exceptions are ZnSe and CuBr, which have PRs that are close to zero. On the other hand, except for CrSe under antiferromagnetic (AFM) state, the selenide compounds with d 3 -d 9 transition metal elements exhibit PPRs that are not sensitive to the magnetic configurations.   Fig. 2b). Therefore, a monoclinic relationship between the PR and the metal-metal distance (or the M-X-M bond angle) is not expected, although there is an indirect correlation between them.
Electronic structure origin of NPR To elucidate the interactions between neighboring M cations, we take VSe, CuI, and AgI as examples and explore their electronic structures in different magnetic configurations. As shown in Fig. 3, the M cation and its four surrounding X anions form a quasiregular tetrahedron, and the crystal field generates lower-energy doubly degenerate states with e (d xy , d z 2 ) symmetry and higherenergy triply degenerate states with t 2 (d x 2 Ày 2 , d xz , d yz ) symmetry. Among all these five orbitals, only the d xy orbital lobes point to the neighboring M atoms. The NNN interactions among M atoms are thus dominated by the d xy -d xy interactions.
For the ferromagnetic (FM) configuration of VSe, the α spin channel of the d xy orbital is lower in energy by U (on-site Coulomb repulsion) than the β spin channel of the d xy orbital, as shown in the projected density of states (PDOS) and the inset interaction diagram in Fig. 3a. Only the orbitals with the same spin interact with each other, which results in bonding states at lower energies and antibonding states at higher energies. The calculated partial charge densities shown in the inset of Fig. 3a confirm the formation of bonding and antibonding d xy states. At an electron energy site~−2.5 eV, most of the charge density is located in the region between the two V atoms that are connected by two Se atoms through the edge sharing of two tetrahedra, indicating the bonding d xy -d xy interaction. At an energy~1.5 eV, a d xy* state character is observed on individual V atom but almost no charge density is located in the region between any two neighboring V atoms. Note that the Fermi energy is set at zero in the figure. It is known that bonding states result in attractive interactions while the occupation of antibonding states induces repulsive interactions. As can be seen in the PDOS of V atoms on d xy orbitals, only the bonding states with α spin are fully occupied, whereas all the other d xy states remain empty. This specific occupation of d xy states therefore leads to strong attractive interactions between neighboring V atoms.
For the AFM configuration of VSe, the d xy orbitals of V atom in the A-sublattice (α spin) in the α spin channel are lower in energy by U than those in the β spin channel, whereas the opposite is true for those in the B-sublattice (β spin). The density of states and the interaction diagram can be seen in Fig. 3b. The solid and dashed curves correspond to the PDOS of V atoms on d xy orbital at the A-sublattice and B-sublattice, respectively. The V atom at the A-sublattice mainly contribute to the bonding states for α spin and antibonding states for β spin, whereas the V atoms at the Bsublattice mainly contribute to the antibonding states for α spin and bonding states for β spin. The bonding states are fully occupied in both spin channels, whereas antibonding states are unoccupied. Apparently, the state occupations under FM and AFM configurations are different, which may lead to a variation in the interaction strength between neighboring metal atoms. For instance, it is claimed in previous work that the strength of the orbital interaction between the two spin sites is proportional to the hopping integral t for the FM arrangement, but to t 2 /U for the AFM arrangement 28 Fig. 3c, d, the interaction between two neighboring M atoms is relatively weak, leading to NPR or near-zero PR. This set of d 10 2D MX compounds therefore offers a promising material platform for the experimental realization of auxetic effects.
Poisson's ratios under strain Different from bulk compounds, 2D materials can undergo a large amount of strain without destroying their atomic structures 29,30 . Specifically, the buckled structure of MX compounds is an important geometric origin of material flexibility. For a buckled structure, a tensile strain along x-direction induces the increase of both bond angle ∠M 1 XM 2 and bond length M 1 X (and M 2 X). Taking AgI as an example, the AgI bond length along the direction of strain increases by only 6.6% under a 25% tensile strain, indicating great material flexibility under strain.
Beyond the linear response regime, it is expected that the metal-metal distances are considerably larger and the d xy -d xy interactions between neighboring metal atoms become weaker, which introduce nonlinear effects and varying PRs under different strain condition. For AgBr, AgI, CuCl, CuBr, and CuI, the PRs decrease by as much as 0.25 as a response to a 5-30% uniaxial strain (Fig. 4). Interestingly, a transition from small PPR to large NPR in CuCl is observed. In the extremely large strain regime (>20%), the PRs of CuBr and CuI reach their minima slowly, whereas those of CuCl, AgCl, AgBr, and AgI increase after reaching the minima.
To understand the mechanical behavior of MX halides under large strain, we study the changes in lattice constants, M-X bond lengths, and M-X-M bond angles under large strain as shown in Supplementary Figure 7). The change in lattice constant along ydirection, and thus the PR, is determined by two competitive geometric effects: the change in M 3 -X (M 4 -X) bond length and the change in M 3 -X-M 4 angle ∠M 3 XM 4 . Under a tensile strain along xdirection, the increase in metal-metal distance induces a weaker metal-metal interaction, which endows the buckled 2D systems less positive or more negative PRs. On the other hand, in the extreme case that the ∠MXM is close to 180°under very large tensile strain, the buckled MX system will collapse to a fully planar structure. In this case, there is no further room for any change of bond angle ∠M 3 XM 4 and the lattice constant along y-direction would decrease under a tensile strain along x-direction, leading to a positive PR.
Considering these two boundary conditions, it is therefore expected that the PRs of all the 2D MX halides would reach their own minima and then increase with even larger tensile strain along x-direction. The strain under which PR reaches its minimum is thus correlated with the bond angle. The critical angles when the PR reaches the minimum for Ag-based compounds are estimated to be close to 140°, whereas those for Cu-based compounds are~132°(as shown in the inset of Fig. 4), which supports the above analysis. Note that the M-X interaction and its important role in the deformation mechanism of 2D compounds have been pointed out in the study of cation-anion interactions in the quasi-planar tetra-coordinated bonding environment of auxetic 2D titanium mononitride 31 .

Exfoliation energies and stabilities of MX compounds
To evaluate the feasibility of experimental synthesis of MX monolayer structures, exfoliation energies are computed by using bilayer systems and considering the van der Waals interlayer interactions. The computational details are included in the Supplementary Methods and the results are shown in Supplementary Figure 8. For all the MX compounds we studied, the interlayer interactions (estimated by the exfoliation energies) are rather weak (ranging from 0.2 to 20 meV/Å 2 ) in comparison with those of existing monolayer 2D systems (as large as 30 meV/Å 2 ) 32 , which is encouraging for the experimental realization of monolayer MX systems from their bulk structures.
Next, we perform phonon calculations and ab initio molecular dynamics (MD) simulations to evaluate the dynamical and thermal stability of MX compounds. The results on phonon spectra and MD simulations are shown in Supplementary Figures 9-11. Ten out of the 14 compounds (except NiSe, CuCl, CuI, and AgCl) exhibit no imagery frequency in their phonon spectra, indicating dynamical stability. For those compounds that host NPRs at equilibrium lattice constants or under small strain, CuBr, AgBr, and AgI exhibit promising dynamical stability. For one of the most promising candidate systems with NPRs, AgI, the MD simulations show its good thermal stability at room temperature.
In summary, based on a material discovery hypothesis for NPR, we perform a structure data mining and identify a family of monolayer 2D MX materials with buckled structures within the p4mm plane group. Using first-principles calculations based on density functional theory and a hybrid functional, we propose that a set of compounds in this family, including CrSe, CuCl, CuI, AgCl, AgBr, and AgI, exhibit auxetic mechanical behaviors while others have more-usual and PPRs. We show that the fundamental difference in their mechanical behavior originates from electronic structure and can be interpreted by analyzing the bonding/ antibonding interactions between a subset of d states of transition metal atoms at the next-nearest-neighboring level. When bonding states dominate the metal-metal interaction, PPR occurs. The electron occupation of antibonding states induces a weaker metal-metal interaction, which eventually leads to a NPR. The nonlinear strain response in these MX compounds is systematically studied. A strong variation in PR under large strain is observed and explained by a specific deformation mechanism in the large strain regime. We anticipate this set of material Fig. 4 Poisson's ratio of metal halides under large tensile uniaxial strain along x-direction. The general trend is that the Poisson's ratio becomes more negative with larger strain along x-direction until reaching a minimum. In the extremely large strain regime, the further increase of uniaxial strain leads to an increase in the Poisson's ratio. The estimated critical M 1 -X-M 2 angles at which the Poisson's ratios reach the minima are shown in the inset.
predictions in the 2D MX family will promote experimental synthesis efforts as well as the continued development of nanoscale 2D auxetic materials for mechanical applications and beyond.

Density functional theory calculations
Lattice constant and geometric relaxations of 2D MX compounds are performed using density functional theory within projector-augmented wave (PAW) potentials 33,34 and the screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE06) 35,36 as implemented in the Vienna ab initio simulation package (VASP) code 37,38 . A vacuum slab of 20 Å is selected and a plane-wave basis set with an energy cutoff of 600 eV is used. We use a 10 × 10 × 1 Monkhorst−Pack k-point set to sample the Brillouin zone. The structures studied here are fully relaxed until energy and force are converged to 10 −7 eV and 0.001 eV/Å, respectively. Such a high criterion is required due to the sensitivity of in-plane elastic properties to small structure change and the existence of interactions between next-nearestneighboring atoms. For the calculations of magnetic MX compounds with transition metal elements, primitive cells are used and the monolayer structures are fully optimized with a given magnetic configuration (either FM or AFM).

Mechanical properties calculations
To investigate the mechanical properties of monolayer MX structures, we employ a least-squares method [39][40][41] with the formula E S ¼ a 1 ε 2 x þ a 2 ε 2 y þ a 3 ε x ε y under small strain (<2.0%). The PR is calculated by ν ¼ a 3 =2a 1 . The details of the fitting process are given in the Supplementary Information ( Supplementary Figures 2 and 3). In the large strain regime, we apply a tensile strain along x-direction and relax the lattice constant along y-direction (Fig. 4). The PR v is calculated by v = −Δy/Δx, where Δx (Δy) is the deviation of lattice constant along x-(y-) direction from the equilibrium value.