Numerical investigation of the Poisson's ratio of amorphous graphene sheets

The one‐atom‐thick allotrope of graphite, C4 is called monolayer graphene and was discovered in 2004. It is well known for its fantastic electro‐mechanical properties. While transverse contraction and Poisson's ratio were previously studied only for homogeneous or crystalline two‐dimensional (2D) mono‐ and bilayer network structures under uniaxial tensile stress, we now numerically investigate these material properties for amorphous or non‐crystalline monolayer graphene. Here, we pay special attention to the influence of the network heterogeneity. We find the stress‐strain correlation in good agreement with the literature, and we find auxetic behavior, that is, negative Poisson's ratio, over a wide range within the elastic regime.

F I G U R E 1 Monolayer 2D graphene sample with network heterogeneity ℎ = 0.6.2D, two-dimensional.
The influence of the network heterogeneity on the transversal strain and the PR of non-crystalline or amorphous 2D network structures has been investigated for mono-and bilayer 2D Silica by Stratmann et al. [19,20].Now, we numerically investigate transversal strain and PR with respect to the network heterogeneity of amorphous monolayer graphene.

METHODS
In our in-house tool JuMol [21], we simulate MD uniaxial tensile tests with monolayer graphene samples of five different network heterogeneities, 10 samples of each heterogeneity.The samples were generated using the Monte Carlo bond switching algorithm discussed by Bamer et al. [22].
Beginning with a crystalline sample that consists exclusively of rings of size  = 6, a pair of directly bonded C 4 triangles are randomly chosen and flipped.Hereby, the sizes of the adjacent rings change: two opposite rings decrease their sizes, while the two rings in between them increase their sizes by one.
We check the ring neighborhood statistics for unphysical clusters, meaning groups of directly neighbored small rings ( < 6) or groups of adjacent large rings ( > 6), by means of the empirical law of Aboav [23].Furthermore, we check if the bond switch brings the ring size variance  of the network structure closer to a target ring size variance   .The latter is the scaled ring size variance of an experimentally realized 2D Silica sample, measured by Lichtenstein et al. [24],   : In the following, We will call ℎ the heterogeneity factor and we will refer to the sample of Lichtenstein et al. as the reference sample in which they observed ring sizes of  = 4 .. 10 [24].
Only if both criteria are fulfilled, the bond switch is accepted.The potential energy is then minimized.This minimization step brings the network from a purely geometrically consistent (all atoms are fully coordinated, except those located next to the boundaries of the network) back to a physically meaningful structure.
Figure 1 depicts a network sample with heterogeneity of ℎ = 0.6 after a comparatively low number of accepted bond switches.Many of the rings still maintained their original ring size of six shown in grey, while only a few changed their sizes to mostly five (green) or seven (blue), even fewer either decreased their sizes to four (green) or increased them to ring sizes beyond seven (magenta, cyan, and red).Note that larger rings have mostly small neighbored rings, as controlled by the law of Aboav [23].
Ring size distributions of the sample shown in Figure 1, and of the reference sample given by Lichtenstein et al. [24].
Figure 2 depicts the corresponding ring size distribution.Each bar has the same color corresponding to the coloring of rings in Figure 1.The black line is the logarithmic envelope of all bars, as proposed by Büchner et al. [25]: Here,   is a continuous value range around the discrete ring sizes given in Figure 2, for example,   = [3.5, 3.6 .. 10.4, 10.5] which is needed for a continuous envelope for the discrete ring size distribution.The value   represents the mean ring size measured by Lichtenstein et al. [24] which is equal to six.The value  is the ring size variance of the actual ring size distribution to be controlled.Since the depicted network sample already reached its target ring size variance, it is   = 0.6  here.The wider and blurred bars in the background represent the ring size distribution of the reference sample, and the thick, gray line is their logarithmic envelope as given by Equation (2).The smaller ring size variance of the network sample shown in Figure 1, shows up in a much higher, and, much thinner peak of the black graph compared to the gray one.
If we perform many more acceptable bond switches on the network sample of ℎ = 0.6 in Figure 1, we obtain the structure shown in Figure 3 with ℎ = 1.4.Figure 4 now tells us that more than 70% of all rings have changed their size, and the network has a larger number of smaller ( < 6) and larger ( > 6) rings than the reference sample.The corresponding peak of the logarithmic envelope (black) thus is slightly lower, and wider than that of the reference distribution (grey).
During the Monte Carlo bond switching algorithm, we use a pure 2D representation of the network structure.Before starting the actual deformation process, we introduce a small noise into the out-of-plane (-) coordinate of every atom and relax the structure.Figure 5 shows that it forms bumps in -direction, as they can be found in monolayer graphene.
All samples are extended with a strain rate of ε = 0.5 ⋅ 10 −3 over 750 steps according to the athermal quasistatic (AQS) deformation protocol [22].This protocol is defined by three consecutive steps: firstly, a box that contains all atoms is deformed affinely.In our case, the box is stretched in -direction with a given -strain rate while the -length is kept fixed.In the second step, all atoms are displaced affinely, that is, according to their relative positions within the box.In our case, all atoms are only displaced in -direction.In the third step, we minimize the potential energy of the network structure by adjusting the positions of all atoms in three dimensions with respect to their local field of potential energy.This step is necessary to transfer the artificially deformed structure back into its local basin of the potential energy landscape, that is, the non-affine displacements.

F I G U R E 4
Ring size distributions of the sample shown in Figure 3, and of the reference sample given by Lichtenstein et al. [24].
For the third step, an interatomic potential is needed that describes the forces that each atom experiences from its neighbors.We use the Reactive Empirical Bond Order (REBO) potential introduced by Brenner et al. [26]: For each atom , the potential   sums up all repulsion forces of each neighbor atom  in the repulsion potential   and its attraction forces in the attraction potential   , scaled by the bond order term   .The latter takes into account if the connection between both atoms is a single, a double, or even a triple bond.After the -th deformation step, tensile strain in direction  ( = , ) is computed as by using the current and the initial box lengths in direction ,    , and,  0  .We, then, compute the Poisson's ratio as: In the next section, we will consider several quantities □  (□ = , , ) averaged, over all  samples sharing the same heterogeneity.The sample averages will be denoted by a bar overhead, □ .The relation between a quantity □  and its sample average □ is given by

Sample averaged 𝒙-tensile stress σ𝒙 over sample averaged 𝒙-strain ε𝒙
Ebrahem et al. carried out numerical tensile tests to investigate the stress-strain correlation for amorphous graphene monolayers [27].Figure 6 shows our numerically determined average tensile stress σ over average tensile strain ε .The maximum tensile strength σ, is reached by the samples of the lowest heterogeneity ℎ = 0.6 at strain values of about ε ≈ 0.42 and decreases with increasing heterogeneity ℎ.This is in good agreement with the results of Ebrahem et al. [27].At -strain values exceeding 0.35, the graphs start losing their smoothness and begin to oscillate, which can be explained by first energy release events such as bond switches or isolated bond breaks.

Sample averaged 𝒚-strain ε𝒚 over sample averaged 𝒙-strain ε𝒙
In this section, we turn to the relation between average and -strain shown in Figure 7.
During deformation, the sample sets of all heterogeneities show both negative and positive -strain, that is, contraction and expansion.

F I G U R E 6
Tensile -stress σ over tensile -strain ε , both averaged over five sets of each 10 of amorphous monolayer graphene samples with levels of network heterogeneity ℎ = 0.6, 0.8, 1.0, 1.2, 1.4.

F I G U R E 8
Detailed view on the graphs shown in Figure 7 within an average -tensile strain range of 0.0 ≤ ε ≤ 0.02.
A very short period of transverse contraction in -direction ( ε < 0) shows up within the range 0.0 < ε < 0.02.We will have a closer look at the minimum transversal strains within this tensile strain range in Figure 8.
Here and in the following, we will observe that the set of samples with ℎ = 1.0 has an extreme value that becomes smaller or greater the farther the heterogeneity moves away to extreme values (ℎ = 0.6, 1.4).Thus, we can state that network samples with ℎ = 1.0 are in a saturated state.Since they have the same ring size distribution as the reference sample realized by Lichtenstein et al. [24], one may refer to them as benchmark samples in the following and call the effect described above "benchmark sample saturation".
Furthermore, Figure 7 shows a difference between the graphs regarding the lower (ℎ < 1.0) and the larger heterogeneities (ℎ ≥ 1.0).Whereas the three --strain-graphs of the larger heterogeneities form a rather symmetric bow between both zero crossings and reach maximum transversal strains ε, of around 0.15, the graphs of ℎ = 0.6, 0.8 reach their maxima earlier.For further increasing -strains up to ε ∼ 0.2, the graphs decrease with low curvature, forming a downwards-sloping plateau.For higher tensile strains, both slopes approach those of ℎ ≥ 1.0.

Sample averaged Poisson's ratio ν𝒙𝒚 over sample averaged 𝒙 strain ε𝒙
Figure 9 shows the sample averaged PR in -direction for tensile deformation in -direction, ν  , over the sample averaged tensile strain in -direction within the range 0.0 ≤ ε ≤ 0.35.The light gray graph depicts the progress of the sample averaged PR, here exemplary visualized for ℎ = 0.6.As all graphs oscillate heavily, we apply a Savitzky-Golay filter (sgf) with a window length of  = 251 and a polynomial degree of  = 3 to visualize the PR's dependency of the network heterogeneity.So, the five colored graphs in Figure 9 show the filtered, sample averaged PRs, denoted by ν(251,3)  , instead of ν .
After having passed a minimum value below zero near ε ≈ 0.05, the PRs of all five heterogeneities raise again to values greater than zero.Since they are the derivatives of the transversal strains with respect to the tensile strain, the difference between the lower (ℎ < 1.0) and the higher network heterogeneities (ℎ ≥ 1.0) shows up here again: the graphs for ℎ = 0.6, 0.8 make an S-shape and cross those of the higher heterogeneities.This phenomenon is especially pronounced for ℎ = 0.6.Instead, the three graphs for ℎ ≥ 1.0 raise by a slightly raising slope It is important to note that all five PRs exceed values of 0.5.This is no contradiction to classical mechanics because, for 2D materials, the upper limit of PR is 1.0 instead of 0.5, so ν,2, = 1.0 [28].
Figure 10 offers a closer look at the filtered, sample averaged PR within the tensile strain range of the auxetic behavior (0.025 ≤ ε ≤ 0.225).All samples of all heterogeneities show NPRs because they expand in -direction while being F I G U R E 1 0 Detailed view on the graphs shown in Figure 9 within an average -tensile strain range of 0.025 ≤ ε ≤ 0.225.stretched in -direction.Considering the corresponding minimum values, we observe that the absolute minimum value is reached by the samples of ℎ = 1.2 with ν, = ν,,1.2= −0.156at a tensile strain of ε, ν,,1.2= 0.056.
Nevertheless, the minimum value for the benchmark samples (ℎ = 1.0), is not much higher with ν,,1.0= −0.154,occurring at a tensile strain of ε, ν,,1.0= 0.053.If we now assume that the slightly lower minimum PR of the samples with ℎ = 1.2 compared to the one of ℎ = 1.0 is due to some statistical fluctuations by taking into account the relatively low number of 10 samples per heterogeneity set, we can state that the minimum PR increases with ℎ moving away from the value of the benchmark samples, ℎ = 1.0, to either smaller (ℎ → 0.6) or greater values (ℎ → 1.4).And we conclude that the benchmark sample saturation shows up here again, highlighted by the blue arrow in Figure 10.
The zero crossings following the minimum PRs show this phenomenon at first glance: with network heterogeneity moving away from ℎ = 1.0, the zero crossing of each graph shifts to tensile strains, the absolute maximum at ε,0, = ε,0,1.0= 0.165, and the minimum ε,0, = ε,0,0.6= 0.090.We visualize phenomenon by the green in Figure 10.

CONCLUSION AND OUTLOOK
As a groundwork for multiscale membrane modeling of electronic devices, we investigate the influence of the heterogeneity ℎ on the transversal strain and PR of monolayer graphene samples under uniaxial tensile stress.We, therefore, perform molecular dynamics tensile test simulations with five groups of different network heterogeneities ℎ = 0.6, 0.8, 1.0, 1.2, 1.4, 10 samples in each group.
The stress-strain response that we obtain for each heterogeneity is in good agreement with the literature.All five sample sets reveal auxetic behavior, meaning positive transverse strains due to tensile loading, and, as a consequence, NPR.For the correlations of transverse strain and PR over tensile strain, the greatest or lowest extreme values are reached by the heterogeneity of ℎ = 1.0 instead of the minimum (resp.maximum) one and shift to moderate values with ℎ moving to either ℎ = 0.6 or ℎ = 1.4.Since the network samples with ℎ = 1.0 have the same ring size variance as an experimentally realized reference sample, we refer to these samples as benchmark samples and we name the phenomenon as benchmark sample saturation.
In all three correlations (stress, transversal strain, and PR over tensile strain), we observe that the graphs for the lower heterogeneities, namely ℎ = 0.6, 0.8, are different than those of the higher heterogeneities, ℎ ≥ 1.0.This difference is worth being investigated further, together with the benchmark sample saturation.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

F
I G U R E 3 Monolayer 2D graphene sample with network heterogeneity ℎ = 1.4.2D, two-dimensional.

F I G U R E 5
Initial state of a monolayer graphene sample with a network heterogeneity of ℎ = 1.0 before start of the AQS deformation process.AQS, athermal quasistatic.