The normal-auxeticity mechanical phase transition in graphene

When a solid object is stretched, in general, it shrinks transversely. However, the abnormal ones are auxetic, which exhibit lateral expansion, or negative Poisson ratio. While graphene is a paradigm 2D material, surprisingly, graphene converts from normal to auxetic at certain strains. Here, we show via molecular dynamics simulations that the normal-auxeticity mechanical phase transition only occurs in uniaxial tension along the armchair direction or the nearest neighbor direction. Such a characteristic persists at temperatures up to 2400 K. Besides monolayer, bilayer and multi-layer graphene also possess such a normal-auxeticity transition. This unique property could extend the applications of graphene to new horizons.

have not yet been reported. To this end, we herein investigated the NA transition in multilayer graphene and the finite temper ature effects. Unlike the previous study that used REBO potential [7], we employed the AIREBO [15] potential with modified C-C bond cutoffs, which could avoid the unphysical description of mechanical properties at large strains [16][17][18][19][20][21][22]. We firstly reproduced the NA transition in monolayer graphene and further found that the NA transition could remain up to temperature of 2400 K. We then proceeded to examine multilayer graphene and found the similar NA transition behavior. Compared to monolayer, multilayer graphene showed higher ν in the positive regime, while nearly the same value in the negative regime.
We carried out molecular dynamics (MD) simulations for monolayer and multilayer graphene using LAMMPS [23]. The structure of monolayer graphene (figure 1) was configured with the armchair direction as the x direction and the zigzag direction as the y direction. A rectangular 4-atom unit cell ( × 0.426 nm 0.246 nm ) was used to build up all the test samples for convenience. The square monolayer graphene consisted of 35 712 atoms in total with the size of about 300 Å. Such a large sample was prepared in aiming to reduce the fluctuation of system variables. This monolayer sample also served as the building block for bilayer and multilayer graphene with interlayer distance pre-configured as 3.34 Å, which is the same as that in graphite. Two different stacking sequences (AA stacking and AB stacking) were considered for multilayer graphene samples (See table S1 in supplementary information (stacks. iop.org/TDM/4/021020/mmedia)). For each simulation, we first run the isothermal-isobaric ensemble 2 (NPT) dynamics at specified temperatures and zero external pressure for 20 ps to reach equilibrium. The fully relaxed sample was subsequently subjected to uniaxial tensile tests in both the armchair and zigzag directions with a constant engineering strain rate of 10 9 − s 1 in the NVT ensemble. The lattice papameters of fully relaxed samples at 300 K are summarized in table S1 in supplementary information. For all the simulations, Newton's equations of motion were numerically integrated using the velocity-Verlet [24] algorithm with a time step of 0.0005 ps. The temperature and pressure were controlled using Nose-Hoover [25,26] thermostat and barostat, respectively. Periodic boundary conditions (PBC) were applied in both the x and y directions, while a fixed boundary was applied in the z direction. OVITO [27] visualization software was used to generate simulation snapshots.
The original version of the REBO-type potential uses a switching function to cut off the C-C interaction between r cc min (1.7 Å) and r cc max (2.0 Å). Mathematically, the switching function is formulated as, The two parameters take into account the nature of carbon covalent bonding and work well for most  (a) Engineering strain ε y versus ε x for monolayer graphene subjected to uniaxial tensile tests in both armchair (x) and zigzag ( y ) directions at 300 K, respectively. (b) The accompanying ν for tensile tests in the armchair direction, and the NA transition is exclusively observed at about 6% strain. equilibrium structures. However, it is problematic when the inter-atomic distance falls within the switching region under large strains, resulting in artificial strengthening of C-C bonds. To circumvent this problem, many researchers have attempted to adjust the cutoffs. For instance, Sammalkorpi et al [18] [22] validated the usage of 1.92 Å as r cc min in studying pentagon-heptagon defect effects on the strength of graphene. It seems that the aforementioned cutoff modifications are quite system dependent and arbitrary, and therefore lack universal guidance for other researchers. Perriot and coworkers [28] later on proposed a screened environment-dependent REBO potential to replace the explicit switching function with a simple yet efficient screening function to solve this problem. However, as of now it has not been widely acknowledged and incorporated into LAMMPS. It is, therefore, still worthwhile to carry out simulations to test various C-C bond cutoffs for our system. Given the fact that REBO and AIREBO potentials essentially make no difference in simulating monolayer graphene, here we chose AIREBO for the consideration of simulating multilayer graphene in which the long-range interactions among graphene layers play an important role. Based on our test results (See figure S1 in supplementary information), we will use the cutoff of 1.92/2.0 Å for the following simulations. Since the stress-strain curves are extremely insensitive to C-C bond cutoffs at strains less than 10%, we will only cautiously report the results at this strain range. Figure 2 shows the correlation between engineering strain ε y and ε x up to 15% (if applicable) for monolayer graphene subjected to uniaxial tensile tests in both the armchair and zigzag directions at 300 K, respectively. Note that the engineering strain ε y by definition is calculated as ε = − L L L y y y y 0 0 ( ) / , with L y0 and L y as the initial and deformed sample lengths in the y direction, respectively. The black curve is the reproduced data of [7] for comparison. It clearly shows that the red curve overlaps the black one up to 10% strain in the armchair direction. An obvious NA transition occurs at about 6% strain as reported. Interestingly, this NA transition behavior is not observed in the zigzag direction as indicated by the blue curve. Figure 2(b) shows the accompanying ν evolution for tensile tests in the armchair direction according to the classical definition ν ε ε = − y x / . The data is obtained by firstly carrying out the 4th polynomial regression with the red curve shown in figure 2(a), and then taking the first derivative of the fitted function. It shows that ν decreases from around 0.21 to −0.03 when strain ε x increases from 0% to 10%. In comparison with the data reported in [7], a noticeable discrepancy exists at equilibrium (zero strain) in which ν was reported to be about 0.31. Considering the first principle calculation of monolayer graphene in our previous work [29] in which ν was calculated to be about 0.18, the authors might have significantly overestimated ν at the initial stage of the tensile test.
We also investigated the temperature effect on the NA transition behavior for monolayer graphene. Figure 3(a) shows the correlation between ε y and ε x at the temperature range from 50 K to 2400 K. In general, ε y decreases with the increasing temperature, and the occurrence of the NA transition remains at nearly 6% strain in the armchair direction until it goes up to 2400 K, at which ε y remains nearly constant when ε x goes beyond 6%. The same trend is also clearly demonstrated in figure 3(b) that shows the corresponding ν for samples at different temperatures. Interestingly, ν decreases with increasing temperature in the PPR regime, while exhibits completely opposite trend in the NPR regime. In particular, monolayer graphene no longer has the NPR regime at 2400 K.
We examined the effect of thermal strains on the NA transition behavior for monolayer graphene. The lattice constant of graphene which is the distance of a carbon atom to its second nearest neighbors are shown in the figure S2 and detailed in table S2 in the supplementary information as a function of the temperature ranged from 5 K to 2400 K. We noticed a thermal contraction of the lattice constants up to 900 K without external applied mechanical stresses. The total amount of such a thermal contraction is 0.0048 Å. The maximum thermal strain is 0.2%, which might contribute 3.3% to the NA transition. This trivial contribution from thermal strain implies the NA transition is insensitive to the temperature, which is verified in figure 3.
In order to further investigate whether this special behavior is exclusive to monolayer graphene, we extended the study to multilayer graphene, including bilayer, trilayer, tetralayer, and even graphite. For simplicity, A and B denote the stacking element (monolayer graphene) for building up multilayer graphene. Specifically, AA stacking and AB stacking were considered for bilayer graphene; AAA stacking and ABA stacking were considered for trilayer graphene; AAAA stacking and ABAB stacking were considered for tetralayer graphene; and bulk graphite was simulated as six AB stacking bilayers of graphene with PBC. Similar to figure 2(a), figure 4(a) shows the correlation between engineering strains ε y and ε x for samples subjected to uniaxial tensile tests in the armchair direction at 300 K. There are a series of minima of transverse strain ε y when the longitudinal strain ε x increases in all the examined systems, with different values but at the same position ε = 0.06 x . This indicates that the NA transition also occurs at about 6% strain for multilayer graphene. However, the transverse strain ε y decreases much faster with respect to an increase of ε x in multilayer graphene as compared to monolayer graphene in the PPR regime, suggesting that multilayer graphene has larger ν. We also found that the change of the transverse strain ε y increases with almost the same pace in the NPR regime. Moreover, the maximum lateral contraction (ε y ) increases slightly from bilayer AB stacking, to trilayer ABA stacking, to tetralayer ABAB stacking, and eventually saturated for graphite. Interestingly, the opposite trend is observed for AA type of stacking, ε y decreases slightly from bilayer AA stacking, to trilayer AAA stacking, and further to AAAA stacking. Note that the ε y difference among multilayer graphene is very trivial so that the aforementioned trend might not exist at all. However, the difference between monolayer one and multilayer ones is dramatically different.
Similarly, figure 4(b) shows the corresponding ν for the above systems. It can be seen that multilayer graphene has slightly higher ν as compared to that of monolayer one when ε x is below 0.03. The ν of AB stacking graphene (AB bilayer and ABAB tetralayer) slightly increases with increasing layers, and falls behind AA stacking (AA bilayer and AAAA tetralayer) at the same strain. When ε x is above 0.03, all of the systems have practically identical ν, and afterwards gradually slip into the NPR regime at around 6% of ε x . Since 6% of strain is still within the range of elastic domain [30], such a NA transition is reversible upon unloading. As a result, there are great potentials in designing and tuning the graphene-based devices.
The auxeticity is highly desirable in some applications including impact mitigation, sealants, water desalination [31]. Auxetic materials have various applications in stents, skin grafts, smart bandage, artificial blood vessels, batting pads, smart sensors, and aero engine fan blades [32]. Combined with high strength and biocompatibility [33], the unique NA mechanical phase transition could extend the applications of graphene to new horizons.