Temperature-dependent negative Poisson’s ratio of monolayer graphene: Prediction from molecular dynamics simulations

1 Introduction

Graphene has attracted significant interest from scientific community and industry due to its extraordinary physical and chemistry properties [1] since it was first discovered in 2004 [2]. Moreover, scientists have observed the fractional quantum Hall effect [3] and unconventional superconducting phenomenon [4] in graphene. As one of the strongest natural materials ever found, the potential mechanical applications of graphene are also of very high expectations [5, 6, 7, 8, 9]. The thermal contraction (negative thermal expansion coefficients, abbreviated NTECs) has been both experimentally observed [10, 11] and theoretically predicted [12, 13, 14]. However, for the negative Poisson’s ratios (NPRs) of monolayer graphene, there exists disagreement in academia due to conflicting results reported in the open literature.

In general, the Poisson’s ratio is used to describe the phenomenon in a material whose deformation perpendicular to the loading direction. Assuming that a load (or strain) is applied along the 11 direction and a strain along the 22 direction (perpendicular to the 11 direction) will be observed. Hence, the definitions of the total and incremental Poisson’s ratios respectively can be defined as

where strain ε is defined as a physical quantity describing the ratio of deformation to the original length. The NPR is also known as auxeticity [15] which is a highly desirable property for applications in the current and future industry. An early research [16] in 2009 performed a Monte-Carlo (MC) simulation which revealed that graphene may possess NPR when the temperature is over approximately 1700K. On the other hand, researchers have reported that when the applied strain reaches certain high level, the intrinsic auxeticity occurs in the monolayer graphene sheet based on their studies using the molecular dynamics and statics simulations [17, 18] and the first principle theory [19]. It is worth noting that the NPRs observed in these researches [17, 18, 19] are only the incremental Poisson’s ratios, in other words, the total Poisson’s ratios of graphene in their studies are still positive. From the view of classical molecular structural mechanics [20, 21, 22] where a C-C bond is modelled as a beam, pristine graphene may have the total NPR in some special cases depending on the parameters for the bond model. Results from the membrane theory [23, 24] showed that both the total and incremental NPRs can be achieved as large as −1 and −1/3, respectively, when the size of graphene is large enough and the tensile load is low enough.

The MD simulation has certain advantages that simulates an ideal object without any defect and impurity and creates an ideal environment where electricity, thermology and chemistry are isolated which the experimental method may not be able to achieve. However, discrepancies exist in the MD studies on the Poisson’s ratio of suspended graphene due to the selection of different potentials for the C-C bond. There are two main semi-empirical potentials for carbon atoms: reactive empirical bond order (REBO) [25] and Tersoff [26]. Based on these two potentials, many extended potentials such as REBO-II [27], AIREBO [28], REBO2-S [29], Extended Tersoff [30] and Tersoff-S [31], have been developed over the last two decades. From existing literature [32, 33], the total NPR (−0.158) can be obtained by MD simulations with Tersoff potential, but cannot be achieved with the potentials based on REBO. However, the Poisson’s ratio of monolayer graphene cannot be directly measured in experiments and the existing indirect measured results from experiments [34, 35] are of large discrepancies varying from 0.12 to 0.36, which may be caused by using different substrates. Hence, researchers still have doubts on the existence of NPR for monolayer graphene and this motivates the current investigation. In the present study, we will carry out MD simulations using the LAMMPS package [36] to study the mechanical behaviors of monolayer graphene in thermal environments and the Erhart-Albe hybrid potential [37] is selected as the force field of carbon atoms. It is worth noting that the Erhart-Albe hybrid potential [37] between a pair of carbon atoms is built based on REBO [25] and Tersoff [26] potentials.

2 Materials and methods

The MD simulation in modelling a monolayer graphene sheet is illustrated in Figure 1. As depicted in Figure 1(b), the rectangular cell of a monolayer graphene consists of 4 carbon atoms. The distance d₀ between a pair of nearest carbon atoms is set as 1.42Å. Accordingly, the length l_a and width l_b of the cell can be obtained. To avoid the scale effect, the in-plane boundary condition is set to be periodical. Hence, the distance between the margin atoms and simulation box boundary should be strictly controlled to construct perfect hexagons of carbon atoms with the periodical boundary, as shown in Figure 1(c). It is worth noting that 11 and 22 are used to denote the zig-zag and the armchair directions respectively in this study.

Figure 1

MD modelling method of a sheet of free standing monolayer graphene and schematic illustrations of (a) graphene; (b) cell geometry; (c) boundary control

The Erhart-Albe hybrid empirical bond-order potential [37] is employed here. The cohesive energy of this potential is written as a sum over individual bond energies, where V_R and V_Λ are the pairwise attractive and repulsive contributions, respectively, f_C is the cutoff function, b is the bond-order and r is the distance between two carbon atoms. The detailed definitions of these physical quantities can be found in the literature [37].

The MD simulation is implemented in an isothermal-isobaric environment with a Nose-Hoover thermostat which provides good conversation with energy and leads to less fluctuation in temperature. The time step is set as 1 fs. Before the phase of tensile loading, the relaxation is performed by at least 40000 steps to stabilize the temperature for the whole system. For the computation of the TEC, the standard relationship is always used

where α is the TEC with unit of K⁻¹. T and L are the simulation temperature and the corresponding length of the graphene, respectively. T₀ and L₀ are the reference temperature and the length, respectively. The applied strain up to 10% is controlled in the simulation of tension.

3 Results and discussion

The MD simulation of a monolayer graphene with 200.26Å × 200.21Å is performed at a temperature range from 0.1K to 1000K. For the purpose to validate the present MD method, the NTECs with the increment of temperature are calculated and are compared with the existing experimental and numerical results [10, 11, 12, 13, 14] in Figure 2. It should be clarified that the TECs in zig-zag and armchair directions (defined as 11 and 22 respectively) are very close to each other but not the same. It can be seen in Figure 2 that the present NTEC agrees well with the one observed in the experiment [11] below 200K. From Figure 2, it can be found that NTEC values can reach to −1.48×10⁻⁵/K. It is worth noting that the TECs dramatically decline in this range and meanwhile wrinkles on graphene surface occur from absolute zero. From 200K to 1000K, most TEC results shown in Figure 2 have an increasing trend. Hence, we believe that graphene is a thermal expansion material but the wrinkles on graphene surface cause the NTEC in a range of temperature. As demonstrated by the trends shown in Figure 2, the TEC will be of a positive value at a fairly high temperature.

Figure 2

Comparison of TECs for perfect monolayer graphene. The solid line and circles represent α₁₁ and α₂₂ obtained from present method; “dash dot dot” and “short dash” lines correspond to existing MD results [12, 13, 14]; the “dash and dot” lines are from the experimental measurements [10, 11]

The geometry of a monolayer graphene is the same as the one used in the study of TEC and tensile simulations of the graphene in the zig-zag direction and armchair direction are performed at 0.1 K, 100K, 300K and 500K. Figures 3 and 4 show the total and incremental Poisson’s ratios in the zig-zag direction and the armchair direction respectively at 0.1 K, and Figure 3(a) and Figure 4(a) show the detailed resultant strain versus the applied strain in both the zig-zag direction and armchair direction. When the temperature is very close to absolute zero, the carbon atoms in the whole system almost do not vibrate, which can prevent the simulation error caused by particles’ self-vibration. Furthermore, there is no wrinkle on graphene surface at this temperature and then graphene can be treated as a real two-dimensional material. In this condition, the NPRs no matter for total or incremental ones can be observed in both in-plane directions. It is also observed that by increasing the applied strain in one direction, both ν₁₂ and v₂₁ are decreased and the more significant auxetic behavior exhibits in the zig-zag direction. The gap between ν₁₂ and v₂₁ can be continuously increased as large as 27.3% once the 10% deformation is achieved. This phenomenon also provides a proof from another aspect that the graphene sheet is mechanically anisotropic. Still, we observe that the total and the incremental NPRs in the armchair direction exhibit nonlinear variations as the applied strain increases.

Figure 3

NPR of monolayer graphene at 0.1 K when strain is applied along the zig-zag direction. (a) The resultant strain ε₂₂ versus the applied strain ε₁₁. (b) The total and incremental Poisson’s ratios obtained from Figure 3(a)

Figure 4

NPR of monolayer graphene at 0.1 K when strain is applied along the armchair direction. (a) The resultant strain ε₁₁ versus the applied strain ε₂₂. (b) The total and incremental Poisson’s ratios obtained from Figure 4(a)

To investigate the thermal effect on the NPR in graphene, the resultant strain and the total NPR curves during the tension process at various temperatures of 100 K, 300 K and 500 K are respectively illustrated in Figures 5 and 6. The applied strain for the case in Figure 5 is in the zig-zag direction while the applied strain is along the arm-chair direction in Figure 6. It is observed that the amplitudes of fluctuations caused by the system self-vibration at 300 K and 500K are very close to each other and they have larger fluctuations than those at 100 K as shown in Figure 5(a) and Figure 6(a). It is interesting to note that the values of NPRs at 300 K and 500 K are also close to each other when the applied strain is over 6%. As the temperature increases, the ν₁₂ is still larger than ν₂₁, which is also observed from the case at 0.1 K. It means that the NPR is indeed an intrinsic property of graphene and this property is further enhanced by the existence of graphene wrinkles. Unlike other MD results [32, 33], the NPR in graphene is not a constant and varies as temperature and applied strain change in our study. Table 1 shows the detailed values of

Figure 5

The tensile behaviors of graphene at 100 K, 300 K and 500 K in the loading direction of zig-zag. (a) The resultant strain ε₂₂ versus the applied strain ε₁₁. (b) The total Poisson’s ratios under various temperatures

Figure 6

The tensile behaviors of graphene at 100 K, 300 K and 500 K in the loading direction of armchair. (a) The resultant strain ε₁₁ versus the applied strain ε₂₂. (b) The total Poisson’s ratios under various temperatures

NPRs in graphene with 10% applied strain It is found that the present NPR results are much lower than that of -0.158 as predicted in the open literature [32, 33].

4 Conclusion

This paper deals with the thermal effect on the NPR and NTEC of monolayer graphene by MD simulation with a hybrid potential. The boundary conditions of the graphene are set as periodical in both the zig-zag and armchair directions. In this study, we find that the TEC maintains negative up to 1000 K and according to the trend the TEC may become positive when the temperature is increased further. The existing MD results indicate that NPR can be only achieved in graphene by using Tersoff potential, but we find that NPR of graphene can also be obtained by using anew one, i.e. Erhart-Albe potential. It is also found that the NPRs including both absolute and differential ones in the armchair direction show strongly nonlinear variations in the tensile process when the temperature is close to absolute zero. Unlike TEC, the auxeticity of graphene becomes more and more remarkable when temperature is increased. Hence, we conclude that NPR is an intrinsic characteristic for graphene.