Discontinuous evolution of the structure of stretching polycrystalline graphene

TL;DR

This paper reveals that polycrystalline graphene exhibits abrupt, avalanche-like out-of-plane displacements and hysteresis under stretching, affecting its mechanical and thermal properties, as shown through computer simulations.

Contribution

It demonstrates the discontinuous, hysteretic response of polycrystalline graphene's structure to stretching, highlighting a new aspect of its dynamical elasticity.

Findings

Discontinuous displacements observed in simulations

Hysteretic behavior during stress cycling

Implications for mechanical and thermal properties

Abstract

Polycrystalline graphene has an inherent tendency to buckle, i.e. develop out-of-plane, three-dimensional structure. A force applied to stretch a piece of polycrystalline graphene influences the out-of-plane structure. Even if the graphene is well-relaxed, this happens in non-linear fashion: occasionally, a tiny increase in stretching force induces a significant displacement, in close analogy to avalanches, which in turn can create vibrations in the surrounding medium. We establish this effect in computer simulations: by continuously changing the strain, we follow the displacements of the carbon atoms that turn out to exhibit a discontinuous evolution. Furthermore, the displacements exhibit a hysteretic behavior upon the change from low to high stress and back. These behaviors open up a new direction in studying dynamical elasticity of polycrystalline quasi-two-dimensional systems, and in particular the implications on their mechanical and thermal properties.

Equations

$$E_0$$

$$E_{\sigma }=E_0-\sigma L_x{L}_y=E_0-\sigma S$$

$$P=\min \left\{1,\exp \left(\frac{E_b-E_a}{k_{B}T}\right)\right\}$$

$$A=\frac{⟨{\left({\mathbf{r}}_i-{\mathbf{r}}_{i,A}\right)}^2⟩}{L_x{L}_y}$$

Full text

Introduction. Graphene, a crystal of carbon atoms arranged in a honeycomb lattice, shows a plethora of exotic mechanical and electronic properties, and emerged as a paradigmatic example of a crystalline membrane embedded in the three-dimensional space Castro Neto et al. (2009); Geim (2009); Nair et al. (2012); Smith et al. (2013); Dolleman et al. (2016); Cartamil-Bueno et al. (2016); Lee et al. (2008); Milowska et al. (2013). In particular, it exhibits an intrinsic tendency to spontaneously buckle when it is polycrystalline, i.e, when it features many crystalline domains, or due to the presence of lattice defects, such as disclinations and Stone-Wales defects Fasolino et al. (2007); Banhart et al. (2011); Zhang et al. (2014a). This effect arises due to the competition of stress, introduced either by mismatch between the domains or by the defects, and the tendency of the membrane to relieve the stress by bending, which yields a rather rich landscape of configurations for a relaxed membrane Witten (2007). In addition, a graphene sheet can experience external stress, for instance when graphene is held by clamps which exert a pulling force, which further enriches the landscape of the ground state configurations.

It is well known how a polycrystalline or defected graphene sheet relaxes when subjected to a constant, static external stress Si et al. (2016), but the case of dynamic strain remained rather unexplored. The latter case can be experimentally relevant as external perturbations creating strain are in reality time-dependent. In particular, when external stress is gradually applied to a polycrystalline graphene sample, it is of a fundamental and practical importance to establish whether the change of the shape of the graphene membrane is continuous or it follows a discontinuous path in this rather complex configuration space.

In this paper, by performing computer simulations, we show that the evolution of the shape of a polycrystalline graphene membrane is discontinuous: as the stress is uniformly increasing, occasionally a tiny increase in stretching force induces a significant displacement, analogous to avalanches, see Fig. 1 and Supplemental Material at [URL will be inserted by publisher] for a video of the evolution of the sample due to the external stress that highlights this avalanche-like behavior, which, as a result, can create vibrations in the surrounding medium. Furthermore, if the stretching force is then decreased again to its starting value, the system does not follow the same path in the configuration space, but rather exhibits a hysteretic behavior while undergoing a cycle in the configuration space, see Fig. 2. The change of the profile of the sample takes place through the creation and annihilation of ridges and vertices, with the elastic energy concentrated in these defects, as shown in Fig. 3. Our findings should have implications for the dynamical elasticity in polycrystalline and defected sheets of graphene, and, in particular, we expect that various elastic moduli, phonon density of states, and thermal conductivity will be affected.

Model. We opted for the empirical graphene potential introduced in Ref. Jain et al. (2015a), which is based on Kirkwood’s potential Kirkwood (1939). This potential has been used, for instance, for studying the long-range relaxation of structural defects Jain et al. (2015a), for probing crystallinity of graphene samples via their vibrational spectrum Jain et al. (2015b), for the study of twisted and buckled bilayer graphene Jain et al. (2017a) as well as of the shape of a graphene nanobubble Jain et al. (2017b). In detail, for a two-dimensional hexagonal network for which out-of-plane deformations are also allowed, this potential can be written as

[TABLE]

with $d=$1.420\text{,}\mathrm{\SIUnitSymbolAngstrom}, the ideal bond length for graphene. A two-body bond stretching energy contribution is parametrized by $\alpha=$26.060\text{\,}\mathrm{eV}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}^{2}, $\beta=$5.511\text{,}\mathrm{eV}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}^{2} controls the bond shearing contribution, while $\gamma=$0.517\text{\,}\mathrm{eV}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}^{2} corresponds to the energy cost for the out-of-plane deformation of the graphene membrane Jain et al. (2015a). Since the periodic boundary conditions apply, we can represent the contribution to the elastic energy due to external stress by an additional term of the form

[TABLE]

where $L_{x}$ and $L_{y}$ are the lateral dimensions of the two-dimensional periodic box, $S$ its surface and $\sigma$ the parameter that controls the strength of the stretching force.

We first generate a polycrystalline flat sample by following the procedure described in Ref. Jain et al. (2015a). A sample is allowed to relax to a lower energy configuration following the Fast Inertial Relaxation Engine (FIRE) algorithm Bitzek et al. (2006). The values of the parameters in this algorithm ($N_{min}$, $f_{inc}$, $f_{dec}$, $\alpha_{start}$ and $f_{\alpha}$) are taken as suggested in Ref. Bitzek et al. (2006). Further relaxation of the sample requires topological changes in the network, which we perform through so called bond transpositions Wooten et al. (1985). In this procedure, four connected atoms are selected, two bonds are then broken and reassigned between them. In the implementation we use an improved algorithm Barkema and Mousseau (2000) that avoids complete relaxation before rejecting a bond transposition. After a random bond switch the complete sample is relaxed and the new configuration is accepted with the Metropolis probability

[TABLE]

which includes the effect of thermal fluctuations: a move that increases the energy of the sample might be accepted, with a higher probability as temperature increases. However, note that no bond transposition is allowed while we perform the cycle in the configuration space. Temperature does therefore partially influence the evolution of the sample, and as long as there are still topological lattice defects around which the stress can accumulate, the discontinuous evolution persists.

Once this flat sample is sufficiently relaxed, every atom is placed at a random non-zero distance from the two-dimensional plane and allowed to relax to a buckled three-dimensional configuration. As topological defects increase the elastic energy of the sample, both globally and locally, they cause stress in the material that, once the material is allowed to relax in a three-dimensional configuration, is released by buckling Jain et al. (2015a). Using such configuration of the graphene membrane, we perform a cycle in the configuration space by manipulating the stretching force $\sigma$ in the empirical potential, given by Eq. (2), while the topology of the sample is kept fixed. During this process, we increase the stress in steps of size $\Delta\sigma$, and, after each step, the sample is allowed to relax. Once the value of the maximum stress $\sigma_{max}$ is reached, it is decreased in the same way. In this work, $\Delta\sigma=$8\text{\times}{10}^{-9}\text{,}\mathrm{u}\text{,}{\mathrm{\SIUnitSymbolMicro s}}^{-2} and $\sigma_{max}=$8\text{\times}{10}^{-6}\text{\,}\mathrm{u}\text{\,}{\mathrm{\SIUnitSymbolMicro s}}^{-2}, with $\mathrm{u}$ the atomic mass unit.

Between each stretching step, the atoms in the graphene membrane translate to a different position, and the relevant translations are only the non-affine (intrinsic) ones, which are defined in the reference frame fixed to the sample itself, and therefore they are not related to the expansion or contraction of the sample. In order to characterize the evolution of the graphene sample under the applied external stress, we thus use the non-affinity parameter defined in Ref. Huisman et al. (2010) as

[TABLE]

where ${\bf r}_{i}$ is the in-plane position of the atom $i$, while ${\bf r}_{i,A}$ is its expected position due to the expansion of the reference frame and $\langle\cdot\rangle$ is the average on all the atoms.

We can also extract the local energy distribution at a given moment with distinct energy contributions for the different terms in the total elastic energy, given by Eq. (1). This is computed by dividing the two body energy of every bond equally between the two atoms and assigning the three body and out of plane energy to the central atom.

Numerical simulations and results. We first adiabatically stretch and relax a sample with $N=1600$ atoms in the way we previously described. After completing a cycle in the parameter space, the sample is in a different configuration, as it can be seen in Figure 1. By repeating this procedure a limited number of times, our sample ends up in a stable configuration. After repeating the cycle, the system reaches again a stable configuration, which is, however, different than the initial one.

The same procedure is then applied to a large sample containing $N=3200$ atoms but in this case we further characterize the sample during a gradual application of the external stress.

We start with a stable configuration and perform a cycle during which we compute the non-affinity parameter, defined in Eq. (4) and its derivative with respect to the stretching parameter. Figure 2 shows that the evolution of the sample during the straining is clearly discontinuous. Significant changes in the configuration are caused by small increases in stretching force, which results in discontinuities in the derivative of the non-affine parameter.

We now analyze in more detail these discontinuities, in particular how the local energy distribution changes. As we can see in Fig. 3, the change of energy is not uniform in the sample. The energy difference, both positive and negative, is concentrated on and around the ridge defects in the sample Witten (2007); Fasolino et al. (2007). Separating the different energy contributions, we notice that, at each transition during the straining, the sample increases its internal energy $E_{0}$ given by Eq. (1), even though not all the terms contributing to it are necessarily positive. Nevertheless, the total elastic energy decreases during these transitions as the contribution of the stretching term, given by Eq.  (2) is negative, indicating a sudden increase of the surface area of the sample. This is confirmed by subtracting the affine component corresponding to the surface expansion or contraction during the transition, computed through a fit of the data points immediately preceding the transition, which show a significant non-affine component (see Supplemental Material at [URL will be inserted by publisher for a breakdown of the affine and non-affine components of the surface changes during the transition). Finally, a minimal size of the sample in which this effect can be observable in our numerical experiment can be estimated as follows. In Fig.  2 the highest peak is $2.7\times 10^{-3}$, while the threshold for the observation of the discontinuity is $0.4\times 10^{-3}$. Given that the sample we used in our simulations (Fig. 3) contains about 15 grains, and assuming a linear scaling of the size of the discontinuity with the number of grains, yields an estimate of a critical size of about 3 crystalline grains.

Discussion & outlook. In summary, we established that the evolution of graphene sheets is discontinuous under gradually applied external stress. It proceeds through a series of avalanche-like processes in which the ridges and vertices are created and annihilated, with the energy concentrated in these defects. Furthermore, the behavior of the graphene membrane is hysteretic: the system does not follow the same path back in the configuration space. Our results imply that if twisting is applied gradually to graphene bilayers and multilayers, the change of the shape of such structures should be discontinuous in nature, which could be possibly relevant to recent experiments Cao et al. (2018a, b); Yankowitz et al. (2019); Lu et al. (2019). Our results should also be pertinent to other two-dimensional polycrystalline materials, including phosphorene, MoS2 and hexagonal boron nitride (h-BN) whose experimental realization has been recently reported in Refs. Pei et al. (2016); Zhang et al. (2014b); Xie et al. (2017), as well as in van der Waals heterostructures Geim and Grigorieva (2013). A continuous application of stress to a graphene elastic membrane necessarily involves a creation and annihilation of the topological defects in the form of ridges and vertices, which can be consequential for the electronic and mechanical properties of these systems. On the other hand, our findings should motivate experimental studies of the dynamical elasticity in these systems, in particular the evolution of the elastic moduli, vibrational density of states and thermal conductivity with the adiabatically applied external stress. Finally, although we considered a specific case of monolayer graphene, we expect that our findings will be applicable to generic membranes embedded in three-dimensional space. In particular, within this approach it would be interesting to study the dynamical properties of origami metamaterials Santangelo (2017).