# Mechanical properties of pristine and nanocrystalline graphene up to ultra-high temperatures

Antonio Gamboa-Suárez, Yenis Seuret-Hernández, Jean-Marc Leyssale

► To cite this version:

Antonio Gamboa-Suárez, Yenis Seuret-Hernández, Jean-Marc Leyssale. Mechanical properties of pristine and nanocrystalline graphene up to ultra-high temperatures. Carbon Trends, 2022, 9, pp.100197. 10.1016/j.cartre.2022.100197. hal-03832930

HAL Id: hal-03832930

https://hal.science/hal-03832930

Submitted on 28 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mechanical properties of pristine and nanocrystalline
graphene up to ultra-high temperatures

Antonio Gamboa-Suárezª, Halis Yenis Seuret-Hernándezª, Jean-Marc Leyssaleᵇ,\*

ªCentro de Investigaciones Químicas-IICBA, Universidad Autónoma del Estado de Morelos, Cuernavaca 62210, Morelos, Mexico

ᵇUniv. Bordeaux, CNRS, Bordeaux INP, ISM, UMR 5255, F-33400 Talence, France

Abstract

Despite being the basic building block of most carbon materials used in
high and ultra-high temperature applications, the mechanical properties of
graphene at such temperatures remain totally unknown. Here we compute
these properties using tensile molecular dynamics simulations at tempera-
tures ranging from room temperature to 5000 K (i. e. about the melting
temperature) for pristine and nanocrystalline graphene. Mechanical simu-
lations at such temperatures are made possible via the use of the screened
environment dependent reactive empirical bond order potential [Perriot et
al. Phys. Rev. B 88, 064101 (2013)]. We show that even at the high-
est temperature, both graphene and nanocrystalline graphene retain high
stiffness, strength and fracture strain. Analyses of tensile curves and tra-
jectories (including an analysis of in-plane non-affine atomic displacements)
show that both pristine and nanocrystalline graphene show brittle rupture,
whatever the temperature, even though some early sign of ductility can be

Email address: jean-marc.leyssale@u-bordeaux.fr (Jean-Marc Leyssale)

<footnote>*Corresponding author</footnote>
detected for the nanocrystalline graphene at 5000 K. This rules out the possi-
bility of a ductile fracture in graphene at any temperature below the melting
point. Finally, the simulations also highlight the importance of bond rotation
(Stone-Wales) based defects in the fracture initiation mechanism for pristine
graphene. Especially at intermediate temperature (3000 K) crack nucleation
in pristine graphene is shown to take place exactly at the moment and place
of the first Stone-Wales defect formation.

Keywords: graphene & polycrystalline graphene, elasticity, failure, ultra
high temperatures, simulation

**1. Introduction**

Owing to their low weight, high melting point and excellent mechan-
ical and thermal properties, carbon-based materials are praised as struc-
tural elements for applications involving high and ultra-high temperatures
[1]. Graphite’s melting point, 4800 ± 100 K according to most recent inves-
tigations [2], is the highest among known materials and much higher than
those of refractory metals (≤ 3700 K) and their carbides (≤ 4300 K)[1].
The room temperature in-plane stiffness (≈ 1 Tpa), strength (≈ 130 GPa)
[3] and thermal conductivity (≈ 5000 W/m/K)[4] of its single layer con-
stituent, graphene, are outstanding and can be combined with the much
weaker out-of-plane values in graphite to design complex architectures as in
carbon fiber reinforced carbon/carbon composites for instance, finding ap-
plications in extreme environments in nuclear energy production as well as
aerospace and aeronautics [5, 6]. Despite their importance, the mechanical
properties of carbon materials at high temperatures are poorly documented
in the literature so far, probably owing to the difficulty of performing such
measurements. Only a few reports have been published in the case of poly-
crystalline graphite [7–10], up to about 3000 K, and carbon fibers [11, 12] up
to about 2300 K. However, experimental values of the mechanical properties
of monocrystalline graphite [13, 14] and graphene [3] are only known at, or
close to, room temperature.

Molecular dynamics (MD) simulation is a particularly well-suited tech-
nique to investigate the mechanical properties of materials. Elastic prop-
erties can be obtained either from the equilibrium fluctuations of strain or
stress during simulations at constant stress or strain, respectively [15, 16], or
from simulations at finite small deformations (compressive, tensile or shear)
[17]. Provided that a suitable interatomic potential is used, allowing to accu-
rately describe bond breaking or formation, non-linear mechanics, including
fracture properties can also be obtained from MD simulations at large de-
formations [18]. Therefore, numerous MD-based predictions of mechanical
properties have been reported for graphene-derived systems.

Grantab et al. [19] first predicted the Young’s modulus (E) and fracture
strength of pristine graphene and graphene with tilt-boundaries by simulat-
ing tensile tests at constant strain rate, followed by many others [20–24]. The
effects of point defects, including Stone-Wales (SW) defects [25–29], single
vacancies (SV) and divacancies (DV) [29, 30] as well as mixtures of differ-
ent defects [31] have also thoroughly been investigated. All these studies
converged in showing that the presence of defects lead to a reduction in frac-
ture strength (σ*) and strain (ε*), and to a lesser extent in the stiffness (E)
of the material, yet the defect concentrations considered remains somehow
unrealistically large, due to the small size of the considered systems (few
nm).

The fracture behavior of several models of polycrystalline graphene, or
nano-crystalline graphene as investigated grain sizes varies from few nm to
few tens of nm, has been reported in numerous papers [18, 32–37]. In all these
studies E was shown to increase and $\epsilon^*$ to decrease with increasing grain size.
The effect of grain size on $\sigma^*$ is less clear as some report a decrease, the pseudo
Hall-Petch effect [18, 34]; others report an increase termed as inverse pseudo
Hall-Petch effect [35, 36]; while a recent report by Han suggest a transition
from pseudo Hall-Petch to inverse pseudo Hall-Petch at grain sizes around 3
nm [37].

Some authors have also investigated the effect of strain rate on the ob-
tained properties [23, 30, 38]. Extrapolating the data to experimental strain
rates (~ 1 s⁻¹), Daly et al. have shown that tensile strength can be reduced
by a up to about 50 %, compared to data at typical MD strain rates (~
10⁹-10⁷ s⁻¹) for graphene with vacancy defects [30].

The effect of temperature on the mechanical properties of graphene has
also been the subject of numerous studies [25, 28, 29, 38–42] and a general
decrease in all properties with increasing temperature was observed. How-
ever, the highest temperature covered was only of 2400 K [38], far from the
melting point of graphene. Common to all these works is the use of the
second-generation reactive empirical bond order (REBO) [43] potential or
its variant including van der Waals interaction, the AIREBO potential [44],
which are intrinsically short-ranged as bonding has to be limited to first
neighbors only. More precisely, bonding interactions are truncated at a cut-
off value $r_c = 2$ Å in these potentials. While a smooth switching function, from 1.7 Å to $r_c$ is normally used in these potentials, the latter has been shown to create artifacts, including an anomalous bump in the stress-strain curve and artificial ductility under tension [19]. A sharp cutoff at, or close to, $r_c$ is thus generally used instead, yet it can lead to artefacts at high temperature and under strain when bond vibrations can reach the cutoff, especially in defective regions.

Modifications of the REBO potential have been proposed to extend the range of the bonding term above 3 Å while ensuring that interactions with second neighbors are not included, through a screening function [45, 46]. Gamboa et al. [18] have used the screened environment dependent REBO potential (SED-REBO) of Perriot et al. [46] to investigate the mechanical properties of graphene and nanocrystalline graphene at room temperature, and compared the properties to those obtained with the REBO potential with a hard cutoff at 2 Å. It was observed that the REBO potential leads to a significant underestimation of both $\sigma^*$ and $\varepsilon^*$, by about 40 %. Building on this work, we present in what follows the mechanical properties of graphene and a polycrystalline graphene model predicted with MD simulations using the SED-REBO potential at various temperatures, from room temperature to 5000 K.

## 2. Methods

Two 2D carbon sheets are considered in this work. The first one is a pristine graphene (PG) sheet of 10032 atoms. The second one, shown in Fig. 1, is a nanocrystalline graphene model (NCG) of similar size (10001 atoms).
<figure><img src="image_2.png" /><figcaption>Figure 1: Top view of the periodic NCG sheet. Rings are color-coded according to the definition found in Ref. [47]: crystalline hexagons (green); non crystalline hexagons (violet); pentagons (red) and heptagons and larger (cyan).</figcaption></figure>
It was prepared [18] through a liquid quench molecular dynamics simulation, using the REBO potential [43], restrained to a plane with an external (cosine-like) potential and further annealed for 25 ns at 4000 K to enhance crystallinity. It contains a large majority of hexagonal rings (about 94 %) arranged in nanometer sized crystallites of arbitrary shapes, separated by both low-angle (discontinuous) and high-angle (continuous) grain boundaries with alternating pentagon/heptagon pairs. The average size of crystallites has been estimated to be between 5 and 8 nm, depending on the calculation method. This is described in Ref. 18, together with all the details on the generation of this model. It is worth mentioning that this model was chosen as it shows highest thermal stability with respect to other models obtained via shorter annealing time. The latter could experience significant structural evolution, mostly an increase in grain size, during high temperature tensile testing. Both systems are enclosed in periodic orthorhombic cells of about 16 nm length in the plane directions and 10 nm in the normal direction.

The mechanical properties of these two carbon sheets are computed at increasing temperatures, from 300 K to 5000 K, using the molecular dynamics methodology described in Ref. 18 and implemented in an “in house” MD code. The interatomic interactions are described with the SED-REBO potential [46], deriving from the second-generation REBO potential [43]. The SED-REBO potential energy is given by

$$E = \frac{1}{2} \sum_{j \neq i} S_{ij} [V_R(r_{ij}) + b_{ij}V_A(r_{ij})] \quad (1)$$

where $V_R$ and $V_A$ are the repulsive and attractive terms, depending only on the distance $r_{ij}$ between atoms i and j; and $b_{ij}$ is an environment dependent bond order term. The most significant difference between SED-REBO and
REBO is the presence of an angular dependent screening function $S_{ij}$ in place of the short-range switching function of the REBO potential. The screening function allows to progressively turn off the interactions when a third atom $k$ get inserted in between $i$ and $j$. Therefore, the interaction range of the potential can be extended much beyond the first neighbor distance. In SED-REBO, $V_R$ and $V_A$ are described using analytical rational functions which, together with their first and second derivatives, smoothly vanish at $r_{ij} = 3.3$ Å. The bond order term, $b_{ij}$, is mostly kept unchanged from the REBO potential, the only difference being the switching function applied to the $b_{ij}^{\sigma-\pi}$ term which is replaced by the product of the switching function $S_{ij}$ by a polynomial $f_c(r)$, ensuring continuous first and second derivatives at the boundaries of the switching domain. A detailed description of the SED-REBO potential is available in Ref. 46. During its implementation, few misprints in the equations and parameters were identified in Ref. 46; these are corrected in the appendix.

Equations of motions are integrated using a velocity Verlet integrator [48] operating with a timestep ranging from 0.5 fs (300 K) to 0.25 fs (5000 K), ensuring a correct energy conservation under equilibrium microcanonical (NVE) runs. The code also uses weak coupling (Berendsen) thermostat and barostat [49]. The barostat allows independent treatment of the diagonal stress (pressure) tensor elements. Time constants of 50 and 250 fs are used for the thermostat and barostat, respectively.

Prior to mechanical testing, the models are relaxed at the suited temperature and zero mechanical stress, with independent control of the stress elements in X and Y directions (the size of the box in the Z direction is held
fixed in all the simulations), for about 50 ps. Stress elements, obtained in
N.m⁻¹ are converting to pressure units assuming a thickness of 3.35 Å for
the graphene sheets. Then, mechanical properties are computed using the
relaxed uniaxial tensile tests (UATT) described in Refs 17, 18, 50. In this
method, a constant engineering strain rate of 10⁹ s⁻¹ is applied to the sys-
tem in one of the in-plane directions while the other is simulated at con-
stant zero stress. Engineering strain (ε_E) is then converted to true strain
ε = ln(1 + ε_E). Therefore, the Young's modulus is directly obtained from
linear regression of the stress-strain curve.

**3. Results and Discussion**

3.1. *Pristine graphene*

Fig. 2 shows examples of the tensile curves obtained for PG along arm-
chair (AC) and zig-zag (ZZ) directions at various temperatures. At 300 K,
the high-strength and high-stiffness with brittle failure behavior already dis-
cussed in Ref. [18] is observed. With increasing temperature, we observe a
progressive decrease in failure strength and strain, as well as a decrease in
modulus. Interestingly, the tensile curves keep a very sharp profile at failure
for all temperatures but 5000 K at which a significantly smoother profile is
observed.

Quantitative values of the mechanical parameters extracted from the ten-
sile curves are given in Fig. 3. It includes values of the Young's modulus
(E), the ultimate strength (σ*) and strain (ε*), defined as the maximum
stress value on the curve and the corresponding strain value, respectively.
Also given is a parameter, γ, characterizing the brittleness of the material
<figure><img src="image_3.png" /><figcaption>Figure 2: Example of tensile curves obtained for PG in (a) AC and (b) ZZ directions at various temperatures.</figcaption></figure>
and defined as $\gamma = \varepsilon_{1/2} - \varepsilon^*$ where $\varepsilon_{1/2}$ is the strain at which the stress has dropped down to half $\sigma^*$ after failure. We note that the definition of $\gamma$ depends both of system size and strain rate, yet allows comparing the different situations encountered in this paper where both size and strain rate are fixed. In the case of an infinitely slow (quasistatic) tensile experiments, or for an infinitely small system, $\gamma \sim 0$ for a perfectly brittle system.

<figure><img src="image_4.png" /><figcaption>Figure 3: Evolutions with temperature of the average (a) Young's modulus E, (b) failure strength $\sigma^*$, (c) failure strain $\varepsilon^*$ and (d) brittleness parameter $\gamma$ (see text for the definition) for PG under AC (squares) and ZZ (circles) tension. Error bars are given to one standard deviation among five independent calculations. Lines are linear fits to the data.</figcaption></figure>

The Young's modulus (E) in the AC direction (Fig. 3(a)) decreases almost perfectly linearly with T in the whole temperature range, from a value of 872 GPa at 300 K down to 591 GPa at 5000 K. Young's modulus in the ZZ direction show slightly lower values, and a less linear dependence on T even though remaining in the same range as AC values. Ultimate strength
($\sigma^*$) decreases perfectly linearly with T from 108 to 39 GPa and 129 to 39 GPa for AC and ZZ directions, respectively. Similarly, the failure strain ($\epsilon^*$) decreases from 0.164 to 0.077 and from 0.238 to 0.078 for AC and ZZ directions respectively.

<table><caption>Table 1: Comparison of the mechanical properties of pristine graphene extrapolated at 0K (linear fits in Fig. 3) with the DFPT calculations of Liu et al. [51]. Note that the armchair and zigzag directions are inverted in Ref. 51; we corrected this in Table 1, for consistency.</caption><thead><tr><th></th><th>DFPT</th><th>This work</th></tr></thead><tbody><tr><td>E<sub>AC</sub> (GPa)</td><td>1050</td><td>889</td></tr><tr><td>E<sub>ZZ</sub> (GPa)</td><td>1050</td><td>806</td></tr><tr><td>σ<sub>AC</sub> (GPa)</td><td>110</td><td>112</td></tr><tr><td>σ<sub>ZZ</sub> (GPa)</td><td>121</td><td>133</td></tr><tr><td>ε<sub>AC</sub></td><td>0.194</td><td>0.164</td></tr><tr><td>ε<sub>ZZ</sub></td><td>0.266</td><td>0.233</td></tr></tbody></table>

The reported values at 300 K, are in fair agreement with the values measured by Lee et al. [3] using AFM: E = 1.0 TPa and $\sigma^* = 130$ GPa. The extrapolated properties at 0 K (Table 1) are also in good agreement with the ideal (0 K) values computed by Liu et al. [51] using density functional perturbation theory (DFPT), especially for the fracture strain and stress. While no experimental data exists so far regarding the mechanical properties of graphene at high temperature, the results in Fig. 3 can still be compared to MD simulations from the literature [25, 29, 38–40, 42]. Although all the available reports indeed show the same trend as observed in this work, namely
a smooth, more or less linear, decrease of fracture strength and strain with
increasing temperature, the rates of these decreases are considerably larger
than the ones reported here. A reduction by 50 % of $\sigma^*$ and $\varepsilon^*$ with respect to
their room temperature values, is obtained at about 2000 K with the sharp-
cutoff AIREBO potential [38, 39] while such a reduction is only obtained
at about 4000 K with SED-REBO, highlighting the importance of properly
accounting for the middle-range (2-3 Å) interactions (note that we compare
to simulations performed using similar strain rates and system sizes to avoid
as much as possible other sources of difference). It is also worth noting that
no simulation data exist in the literature for T > 2400 K, suggesting that
sharp-cutoff potentials could be unstable at high T and under strain.

We notice that all mechanical properties remain at very high values up to
5000 K, despite it is about the melting point for graphene. From the linear
fits in Fig. 3 (b and c) it is tempting to determine bounds for a “mechanical”
melting. Extrapolation of the data in Fig. 3(b) to $\sigma^* = 0$ leads to predicted
mechanical melting points of about 7800 and 7000 K in AC and ZZ directions
respectively. Extrapolation of the data in Fig. 3(c) to $\varepsilon^* = 0$ leads to a similar
estimate in the ZZ direction, 7200 K, and a slightly larger value of about 9200
K in the AC direction. These values are significantly larger than the actual
thermodynamic melting point (see below).

As shown in Fig. 3(d), the $\gamma$ parameter remains close to its 300 K value (~ 0.003) up to about 4000 K, showing only a very slight linear increase. This indicates that the brittle fracture of PG is maintained. However, a significant increase in $\gamma$ is observed at 5000 K (0.01). Nevertheless, even at 5000 K, the material is far from showing the ductile behavior of a metal or
a ductile polymer for which $\gamma$ would be typically larger than one.

<figure><img src="image_5.png" /><figcaption>Figure 4: Snapshots of PG under UATT at 3000 K taken nearby crack nucleation. (a-d): AC direction; (e-h): ZZ direction. Atoms in blue and red in the insets are visual guides to show how crack nucleation is concomitant to SW rotation attempts. Numbers at bottom-left corners indicates the strain level for each frame. Full trajectories are available as supplementary movies M1 and M4.</figcaption></figure>

Visualization of the trajectories at 1000 K and 2000 K did not unravel any difference in the fracture mechanism with respect to earlier reports at room temperature [18]. However, some interesting observations have been made at higher temperatures. Fig. 4 shows snapshots of the crack nucleation at 3000 K taken from two MD trajectories, one AC and one ZZ. In the AC direction we observe that crack nucleation is concomitant to a bond rotation attempt at $\epsilon \sim 0.11$ (see the blue atoms in Fig. 4(a-d)). In the ZZ direction, we observe that a complete Stone-Wales (SW) defect [52] is formed at $\epsilon \sim 0.12$ (Fig. 4(f)) and quickly after ($\epsilon \sim 0.13$) serves as a crack nucleus (Fig. 4(g)).

At 4000 K SW defects start appearing at $\epsilon \sim 0.09$ (Fig. 5). As shown in Fig. 5(a), 4 SW defects, as well as 3 more complex defects based on
<figure><img src="image_6.png" /><figcaption>Figure 5: Snapshots of PG under UATT at 4000 K taken nearby crack nucleation. (a-d): AC direction; (e-h): ZZ direction. Numbers at bottom-left corners indicates the strain level for each frame. Colored atoms indicate the presence of point, mostly SW, defects prior to crack nucleation. Full trajectories are available as supplementary movies M2 and M5.</figcaption></figure>

multiple bond rotations and partially dissociated are observed prior to crack nucleation under AC strain. At $\epsilon \sim 0.098$, two cracks, centered on existing defects, nucleate almost concomitantly (Fig. 5(b,c)) and lead to fracture (Fig. 5(d)). As shown in Fig. 5(e-h), the fracture scenario is very similar under ZZ strain. First, 5 SW defects and an extended defect are formed at $\epsilon \sim 0.094$ (Fig. 5(e)), then a crack nucleates at the extended defect and propagates to complete fracture (Fig. 5(f-h)).

Fig. 6 shows snapshots of PG during UATT at 5000K. As shown in Fig. 6, SW defects appear at much lower strain values of 0.026 (Fig. 6(a)) and 0.039 (Fig. 6(d)) in the AC and ZZ directions, respectively. Then many defects form with increasing strain, up to crack nucleation, which can be located at strains of about 0.075 (Fig. 6(b)) and 0.086 (Fig. 6(e)) for AC
<figure><img src="image_7.png" /><figcaption>Figure 6: Snapshots of PG during UATT at 5000 K. (a-c): AC direction; (d-f): ZZ direction. Numbers at bottom-left corners indicates the strain level for each frame. First SW defects are colored in black in (a,d). Full trajectories are available as supplementary movies M3 and M6.</figcaption></figure>
and ZZ directions, respectively. One note that these crack initiation events
are less marked that what was observed at lower temperatures. As can be
seen in Fig. 6(c,f), the crack profiles show much more disordered edges than
those that were observed at lower temperatures, in agreement with the less
brittle behavior observed from the tensile curves (Fig. 2)

3.2. *Nanocrystalline graphene*

<figure><img src="image_8.png" /><figcaption>Figure 7: Example of tensile curves obtained for NCG in (a) X and (b) Y directions at various temperatures.</figcaption></figure>

Fig. 7 shows examples of tensile curves obtained for the NCG model.
Compared to PG (Fig. 2), we observe a well-documented non-linear elastic
behavior at low strain values [18, 33], associated to flattening of the NCG
which present significant curvature at rest. Besides, an overall reduction of
the tensile properties, compared to graphene, is also clearly visible. Nevertheless, the NCG model also keeps high mechanical stiffness and strength in the considered T range. As for graphene, we also observe a significant reduction in tensile modulus and tensile strength with increasing temperature and a decrease in brittleness at the highest temperatures. However, the evolutions of ultimate strain is less obvious than for graphene. More details are given in Fig. 8.

<figure><img src="image_9.png" /><figcaption>Figure 8: Evolutions with temperature of the average (a) Young's modulus E, (b) failure strength $\sigma^*$, (c) failure strain $\epsilon^*$ and (d) brittleness parameter $\gamma$ for the NCG model under tension along the X (squares) and Y (circles) directions. Error bars are given to one standard deviation among five independent calculations. Lines are linear fits to the data.</figcaption></figure>

Fig. 8(a) shows the evolution of E with T. An almost identical linear decrease is observed in both X and Y directions, from about 800 GPa at 300 K down to about 600 GPa at 4000 K. However, a significant drop in E, to about 400 GPa, is observed when T is increased up to 5000 K, in-
dicating a change in the mechanical behavior of the NCG model. Similar
observations can be made for the ultimate strength σ* (Fig. 8(b)), showing
a linear decrease with increasing T in both directions up to 4000 K (with
slightly more pronounced initial value and slope in the X direction). Then,
as for E, a sudden drop of the value is observed at 5000 K. The behavior of
fracture strain ε* (Fig. 8(c)) shows an almost perfectly linear decrease for
tests performed in the Y direction, from about 8.8 % at 300 K to about 7.7
% at 5000 K. However, the evolution of ε* with T appears more complex in
the X direction. Indeed, whilst ε* values are almost identical to those in the
Y directions for T ≥ 2000 K, values obtained at lower temperatures are sig-
nificantly larger. A possible explanation for this peculiar behavior could be
the extended non-linear elastic behavior observed for the NCG model along
X tension at low temperature (see the tensile curve at 300 K in Fig. 7(a)).
As for PG, the ductility parameter γ shows a slow linear increase with T up
to 3000 K (Fig. 8(d)), albeit with slightly larger values than for graphene,
from about 0.004 at 300 K to about 0.005 at 3000 K. γ then significantly
increases with further increasing T to about 0.025 ± 0.01 and 0.035 ± 0.01 in
X and Y directions, respectively, at 5000 K. This sudden rise in γ indicates
that the NCG approaches the brittle to ductile transition, yet considering
the temperature, melting would certainly occur before such transition could
be observed.

Figure 9 shows the evolution of the NCG model under UATT at 3000
K. At this temperature, the model behaves almost exactly as was described
in Ref. 18 for the same model at room temperature. We observe first the
appearance of many single and few bond cracks, indicated with red arrows
<figure><img src="image_10.png" /><figcaption>Figure 9: Snapshots of the NCG model under UATT (Y axis) at 3000 K. Bonds between defect atoms (i.e. point defects of grain boundaries) in the undeformed NCG model are displayed in black. Numbers at bottom-left indicates the strain level for each frame. Single and multiple bond cracks are highlighted with red arrows and ellipses, respectively. The blue squares are enlargement showing some structural evolutions (Stone-Wales-like) at grain boundaries. The full trajectory is available as supplementary movie M7.</figcaption></figure>
and ellipses, respectively. Single bond cracks are mostly formed by breaking
bonds at grain boundaries that are mostly aligned with the tensile direction
and shared by a hexagon and a heptagon. These defects do not immediately
lead to fracture and can accommodate a significant amount of additional
strain, multiply in number, before larger cracks eventually nucleate (Fig. 9(a-c)). At $\epsilon = 0.088$, an array of small cracks has aligned in a direction almost
perpendicular to the tensile direction (Fig. 9(d)). The latter quickly regroup
into a large unique crack (Fig. 9(e)) that propagates and promote fracture
(Fig. 9(f)). Asides from the crack locations, at 3000 K, the structure of the
NCG is mostly preserved. However, a few bond rotations occurred during
the simulation as shown in the enlargement (blue squares) of two zones in
Fig. 9(a,f). We also observe that most open single and few bond cracks
vanish when away from the main fracture during its growth.

The behavior of the NCG at 4000 K (Fig. 10) is very similar to the one at
3000 K aside from the facts that fracture nucleation and completion occur at
slightly lower strain values. Grain boundary reconstructions are also present
at this temperature, probably even slightly more pronounced than at 3000
K.

At 5000 K, some important differences can be observed (Fig. 11). First,
in addition to the reorganization occurring at grain boundaries, spontaneous
formation and migration of SW defects is also observed in crystalline areas
(green squared area in Fig. 11(a,f))), as was the case for PG at very high
temperatures. Second, a major difference with the 3000 K and 4000 K sce-
narios is that crack propagation occurs on a much longer strain scale: from
the crack nucleation at $\epsilon \approx 0.08$ (Fig. 11(c)) to a state at $\epsilon = 0.15$ (Fig.
<figure><img src="image_11.png" /><figcaption>Figure 10: Same as Fig. 9 for the NCG model under UATT (Y axis) at 4000 K. The full trajectory is available as supplementary movie M8.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Snapshots of the NCG model under UATT (Y axis) at 5000 K. The green square highlights the spontaneous formation (and migration) of a SW defect in a crystalline area. The orange square shows an early stage of melting near the convergence of two fracture fronts threw the periodic boundary (x-axis). The full trajectory is available as supplementary movie M9.</figcaption></figure>
11(f)) at which fracture is not even fully completed. This is a clear sign that we are approaching the brittle-ductile transition. Finally, a last interesting feature observed in this trajectory is the appearance of a melting front at the fracture edge (see the orange squared area in Fig. 11(f)). We note that this melting area seems to be related to the convergence of two fracture fronts through the periodic boundary conditions (on the top and bottom sides of Fig. 11(f)) and is not unique to this system. This will be discussed in more details below.

## 4. Discussion

Local stress or strain analysis close to fracture initiation is a convenient way to characterize materials behavior at fracture. While our current implementation of the SED-REBO do not allow for computing atomic stress elements, we propose in Fig. 12 an analysis of in-plane non-affine displacements of the PG and NCG models close to the fracture initiation. The in-plane non-affine displacement of an atom between is defined as $u_{xy}^{na} = r_{xy} - r_{xy}^{0} - u_{xy}^{a}$ where $r_{xy}$ and $r_{xy}^{0}$ represent the in-plane position of the atom in the current and reference configurations, respectively, and $u_{xy}^{a}$ is the in-plane affine displacement of the atom corresponding to a fully homogeneous (homothetic) deformation. In other terms, $u_{xy}^{na}$ defines the local inhomogeneities of the displacements.

As shown in Fig. 12, no significant difference in the non-affine displacement fields can be observed between PG and NCG at a given temperature. However, some temperature effects are clearly noticeable. At 3000 K, we can see in Fig. 12(a,d) that large $u_{xy}^{na}$ values are mostly observed localized in the
<figure><img src="image_13.png" /><figcaption>Figure 12: Snapshots of the PG and NCG systems around crack nucleation with atoms and bonds color-coded according to the atomic in-plane non-affine displacements. (a-c): PG under UATT in the AC direction at (a) 3000 K, (b) 4000 K and (c) 5000 K; (d-f) NCG under UATT along the Y direction at (d) 3000 K, (e) 4000 K and (f) 5000 K. The reference states for the calculations are the PG at rest and the NCG at $\epsilon = 0.05$ (i.e. in the linear elastic regime obtained after initial flattening of the structure).</figcaption></figure>
neighborhood of the crack nucleus, indicating that strain is homogeneously distributed in the remaining part of the system, as expected for brittle materials. Note that we also observe few hot spots (see Fig. 12(d)) in the NCG around the few non-propagating single of few-bond cracks identified in Fig. 9(c) or at room temperature in Ref. [18]. Increasing the simulation temperature, we observe that $u_{xy}^{na}$ hot spots becomes more extended, indicating a progressive reduction in the brittleness of the materials, in agreement with the increases of the ductility parameter reported above, yet remaining extremely far from values expected for a ductile material ($\gamma \sim 1$).

An important finding of this work is the influence of SW defects on the fracture of PG at high temperatures. Although only two examples (one per tensile direction) per temperature were discussed in the results section (see Fig. 4,5 and 6), we observed the spontaneous formation of SW defects in all performed simulations (5 per temperature and tensile direction) prior to or at the crack nucleation event. The numbers of such defects were of about 1-2 per simulation at 3000 K, 4-6 at 4000 K and slightly more at 5000 K. The spontaneous formation of SW defects, as well as other bond rotation-based defects is extremely unlikely in the absence of strain, even at the temperatures considered here. Indeed, for a PG sheet with $N = 10032$ atoms, the equilibrium number of SW defects can be estimated (neglecting entropic effects) as $N_{SW} = \frac{3N}{2} \exp\left(\frac{E_f}{k_B T}\right)$ where $E_f$ is the defect formation energy and $k_B$ is Boltzmann's constant. Assuming a formation energy of about 5 eV for the SW defect, which is close to both density functional theory calculations [53] and empirical calculations based on REBO-based force-fields [54], one obtains $N_{SW}$ values ranging from about $10^{-80}$ at room temperature to about
$10^{-2}$ at 5000 K. However, it has been shown in prior studies that tensile strain considerably reduces both the SW defect formation energy and the associated energy barrier (~10 eV in the absence of strain) [55]. The presence of one or several SW defects was shown to become energetically favorable with respect to pristine state above a certain strain level [56-58]. Hence, SW defects were considered as plausible precursors to fracture in graphene and carbon nanotubes [58]. Furthermore, performing the simulations at elevated temperatures as those in this work, significantly accelerate barrier crossings. For instance, Leyssale and Vignoles [54] have been able to simulate the diffusion of di-vacancies in graphene at temperatures in the 3000-4000K range, for a process associated to barriers of about 6 eV. The results presented above constitute a clear and definitive proof that SW defects indeed appear under strain, especially at high temperature, and promote fracture in crystalline graphene. SW defects also form in the NCG system, yet without much influence on the fracture properties as crack nucleation occurs at already existing defects or grain boundaries.

A third point to discuss here is the melting domain developing in the NCG system at 5000 K (Fig. 11) during fracture propagation. Even though not presented in section III, focusing on the crack nucleation, partial melting was also observed in other trajectories, both for PG and NCG, for $T \ge 3000$ K. Partial melting can actually be observed for PG under UATT in the AC direction at 3000 K (movie M1), 4000 K (movie M2) and 5000 K (movie M3), in the ZZ direction at 4000 K (movie M5) and 5000 K (movie M6), as well as for UATT for NCG at 5000 K in both Y (movie M9) and X (movie M12) directions. From these movies, it seems that PG seems to be more prone to
melting than NCG, which could result from a stronger shock wave development after fracture in the stiffer PG system. However, a common observation is that the initiation of melting occurs in areas where two misaligned fracture fronts percolate through the periodic boundary, in the direction normal to tension. Care must thus be taken with interpreting this strongly non-equilibrium melting process which could be very sensitive to the simulation conditions, regarding for instance system size, boundary conditions, strain rate, temperature and stress control, etc... Therefore, further work is definitely required to investigate the conditions of fracture-induced melting. In an attempt to locate the melting point $T_m$ of graphene, as described using the SED-REBO potential, 1 ns long MD simulations were performed at various temperatures by steps of 100 K and fixed zero stress. An upper “dynamical” limit of 5600 K was obtained for $T_m$, corresponding to the lowest temperature at which melting was observed. This value is in line with other dynamical estimates based on the LCBOP2 potential [59], about 5000 K [60], or on ab initio molecular dynamics, 4000-6000 K [61]. A more accurate estimate based on nucleation theory for LCBOP2 gives a slightly lower value of 4500 K [59].

Finally, it is important to recall that crack nucleation is a thermally activated process, obeying Eyring's theory, and hence conveys a probabilistic nature [23, 38]. For instance, the characteristic time to fracture nucleation at given stress and temperature reads [23]:

$$ \tau(\sigma, T) = \frac{\tau_0}{n_s} \exp \left( \frac{E_0 - V_a \sigma}{k_B T} \right) \quad (2) $$

where $\tau_0$ is the characteristic time of molecular vibrations in the material, $n_s$ is the number of possible sites for crack nucleation, $E_0$ is the energy barrier to
crack nucleation, $V_a$, an activation volume and $k_B$ is Boltzmann's constant. It immediately follows that this time depends exponentially on $T$. Nevertheless, Eq. 2 shows that $\tau$ also depends on the number of possible crack nucleation sites $n_s$, which scales as the number of atoms in the system, and with the applied stress $\sigma$, which by extension induces a dependence of $\sigma^*$ on strain rate. Here, using the SED-REBO potential we have been able to produce reliable results for the dependence of fracture properties with temperature up to almost melting, which was not possible with shorter ranged interatomic potentials. While system size and strain rate effects on those properties exist in the literature [23, 38], the present work suggest that such investigations should be also performed with an extended range potential as SED-REBO, especially at high temperatures.

## 5. Conclusion

The mechanical properties of pristine and nanocrystalline graphene have been computed using tensile molecular dynamics simulations and the SED-REBO potential. Benefiting from the extended interaction range of the potential, we have been able to characterize both stiffness (Young's modulus) and fracture properties (ultimate strain and stress) up to ultra-high temperature, typically close to the melting point of these materials, which is unprecedented in the literature.

All properties were found to decrease (mostly linearly) with increasing temperature up to 5000 K for PG and 4000 K for NCG. At 5000 K, PG has a Young's modulus of about 600 GPa and ultimate stress and strain of about 40 GPa and 0.08, respectively. NCG also retains impressive properties at
5000 K, with a Young's modulus of about 400 GPa and ultimate stress and strain of about 25 GPa and 0.077, respectively.

Analysis of the tensile trajectories has allowed to finely characterize the fracture mechanism. Brittle rupture was observed for both PG and NCG at all temperatures, with only a few hints of ductility at the highest temperature for the NCG. This work thus rules out the possible existence of a brittle to ductile transition in graphene, as we clearly show here that melting should occur prior to such transition.

The simulations have also shown the importance of the spontaneous formation of SW defects in the crack nucleation mechanism at high temperatures (≥ 3000 K) for PG.

## Appendix

The pairwise repulsive and attractive terms, $V_R$ and $V_A$ are expressed as analytical rational functions:

$$V = \frac{\bar{r}^3}{r^{A_1}} \times \frac{A_2 + A_3\bar{r} + A_4\bar{r}^2 + A_5\bar{r}^3 + A_6\bar{r}^4}{1 + B_3\bar{r} + B_4\bar{r}^2 + B_5\bar{r}^3 + B_6\bar{r}^4} \quad (3)$$

where $\bar{r} = r_c - r$ (instead of $\bar{r} = r - r_c$ in Ref. 46) and $r_c = 3.3$ Å is the SED-REBO cutoff. The parameter $A_1$ was given the wrong sign for both the repulsive and attractive potentials in the original publication. It should read 1.4015028876548725 and 5.8766423956332570 for $V_R$ and $V_A$, respectively. Finally, the switching function $f_c(r)$ should read $1 - \chi^3(6\chi^2 - 15\chi + 10)$ instead of $1 - \chi^3(6\chi^2 - 15\chi - 10)$ for $r_{min} \le \chi \le r_{max}$ where $\chi = (r - r_{min})/(r_{max} - r_{min})$ with $r_{min} = 3.0$ Å and $r_{max} = 3.3$ Å.
[1] Shabalin IL. Ultra-High Temperature Materials I. Springer, Dordrecht; 2014.

[2] Savvatimskiy AI. Measurements of the melting point of graphite and the properties of liquid carbon (a review for 1963–2003). Carbon. 2005;43(6):1115–1142.

[3] Lee C, Wei X, Kysar JW, Hone J. Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene. Science. 2008;321:385–388.

[4] Balandin AA, Ghosh S, Bao W, Calizo I, Teweldebrhan D, Miao F, et al. Superior Thermal Conductivity of Single-Layer Graphene. Nano Lett. 2008;8(3):902–907.

[5] Savage G. Carbon/Carbon composites. London: Chapman & Hall; 1993.

[6] Scarponi C. 13 - Carbon-carbon composites in aerospace engineering. In: Rana S, Fangueiro R, editors. Advanced Composite Materials for Aerospace Engineering. Woodhead Publishing; 2016. p. 385–412.

[7] Malmstrom C, Keen R, Green L. Some Mechanical Properties of Graphite at Elevated Temperatures. J Appl Phys. 1951;22(5):593–600.

[8] Mason I, Knibbs R. Variation with Temperature of Young's Modulus of Polycrystalline Graphite. Nature. 1960;188:33–35.

[9] Smith MC. Longitudinal tensile behavior of h4lm graphite in nitrogen to 2750 °C. J Nucl Mat. 1964;13(1):57–62.
[10] Smith MC. Effects of temperature and strain rate on transverse tensile properties of H4LM graphite tested in helium and in vacuum. Carbon. 1964;1(2):147-153.

[11] Trinquecoste M, Carlier JL, Derré A, Delhaès P, Chadeyron P. High temperature thermal and mechanical properties of high tensile carbon single filaments. Carbon. 1996;34(7):923-929.

[12] Sauder C, Lamon J, Pailler R. Thermomechanical properties of carbon fibres at high temperatures (up to 2000 °C). Compos Sci Technol. 2002;62(4):499-504.

[13] Blakslee OL, Proctor DG, Seldin EJ, Spence GB, Weng T. Elastic Constants of Compression-Annealed Pyrolytic Graphite. J Appl Phys. 1970;41:3373-3382.

[14] Bosak A, Krisch M. Elasticity of single-crystalline graphite: Inelastic X-ray scattering study. Phys Rev B. 2007;75:153408.

[15] Parrinello M, Rahman A. Strain fluctuations and elastic constants. J Chem Phys. 1982;76(5):2662-2666.

[16] Clavier G, Desbiens N, Bourasseau E, Lachet V, Brusselle-Dupend N, Rousseau B. Computation of elastic constants of solids using molecular simulation: comparison of constant volume and constant pressure ensemble methods. Mol Simul. 2017;43(17):1413-1422.

[17] Farbos B, Da Costa JP, Vignoles GL, Leyssale JM. Nanoscale elasticity of highly anisotropic pyrocarbons. Carbon. 2015;94:285-294.
[18] Gamboa A, Farbos B, Aurel P, Vignoles GL, Leyssale JM. Mechanism of strength reduction along the graphenization pathway. Sci Adv. 2015;1(10):e1501009.

[19] Grantab R, Shenoy VB, Ruoff RS. Anomalous Strength Characteristics of Tilt Grain Boundaries in Graphene. Science. 2010;330(6006):946-948.

[20] Wei Y, Wu J, Yin H, Shi X, Yang R, Dresselhaus M. The nature of strength enhancement and weakening by pentagon-heptagon defects in graphene. Nat Mater. 2012;11:759-763.

[21] Rasool HI, Ophus C, Zhang Z, Crommie MF, Yakobson BI, Zettl A. Conserved Atomic Bonding Sequences and Strain Organization of Graphene Grain Boundaries. Nano Lett. 2014;14(12):7057-7063.

[22] Han J, Ryu S, Sohn D, Im S. Mechanical strength characteristics of asymmetric tilt grain boundaries in graphene. Carbon. 2014;68:250-257.

[23] Daly M, Veer Singh C. A kinematic study of energy barriers for crack formation in graphene tilt boundaries. J Appl Phys. 2014;115(22):223513.

[24] Xu N, Guo JG, Cui Z. The influence of tilt grain boundaries on the mechanical properties of bicrystalline graphene nanoribbons. Physica E. 2016;84:168-174.

[25] Wang MC, Yan C, Ma L, Hu N, Chen MW. Effect of defects on fracture strength of graphene sheets. Comput Mat Sci. 2012;54:236-239.
[26] He L, Guo S, Lei J, Sha Z, Liu Z. The effect of Stone–Thrower–Wales defects on mechanical properties of graphene sheets - A molecular dynamics study. Carbon. 2014;75:124–132.

[27] Rajasekaran G, Parashar A. Molecular dynamics study on the mechanical response and failure behaviour of graphene: performance enhancement via 5–7–7–5 defects. RSC Adv. 2016;6:26361–26373.

[28] Juneja A, Rajasekaran G. Anomalous strength characteristics of Stone–Thrower–Wales defects in graphene sheets – a molecular dynamics study. Phys Chem Chem Phys. 2018;20:15203–15215.

[29] Li M, Deng T, Zheng B, Zhang Y, Liao Y, Zhou H. Effect of Defects on the Mechanical and Thermal Properties of Graphene. Nanomaterials. 2019;9:347.

[30] Daly M, Reeve M, Veer Singh C. Effects of topological point reconstructions on the fracture strength and deformation mechanisms of graphene. Comput Mater Sci. 2015;97:172–180.

[31] Xu L, Wei N, Zheng Y. Mechanical properties of highly defective graphene: from brittle rupture to ductile fracture. Nanotechnology. 2013;24:505703.

[32] Kotakoski J, Meyer JC. Mechanical properties of polycrystalline graphene based on a realistic atomistic model. Phys Rev B. 2012;85:195447.

[33] Zhang T, Li X, Gao H. Flaw Insensitive Fracture in Nanocrystalline Graphene. Nano Lett. 2012;12:4605–4610.
[34] Song Z, Artyukhov VL, Yakobson BI, Xu Z. Pseudo Hall-Petch Strength Reduction in Polycrystalline Graphene. Nano Lett. 2013;13:1829-1833.

[35] Sha ZD, Quek SS, Pei QX, Liu ZS, Wang TJ, Shenoy VB, et al. Inverse Pseudo Hall-Petch Relation in Polycrystalline Graphene. Sci Rep. 2014;4:5991.

[36] Mortazavi B, Cuniberti G. Atomistic modeling of mechanical properties of polycrystalline graphene. Nanotechnology. 2014;25:215704.

[37] Han J. The transition from an inverse pseudo Hall-Petch to a pseudo Hall-Petch behavior in nanocrystalline graphene. Carbon. 2020;161:542-549.

[38] Zhao H, Aluru NR. Temperature and strain-rate dependent fracture strength of graphene. J Appl Phys. 2010;108(6):064321.

[39] Zhang H, Duan Z, Zhang X, Liu C, Zhang J, Zhao J. Strength and fracture behavior of graphene grain boundaries: effects of temperature, inflection, and symmetry from molecular dynamics. Phys Chem Chem Phys. 2013;15:11794-11799.

[40] Dewapriya MAN, Rajapakse RKND. Molecular Dynamics Simulations and Continuum Modeling of Temperature and Strain Rate Dependent Fracture Strength of Graphene With Vacancy Defects. J Appl Mech. 2014;81:081010.

[41] Anastasi AA, Ritos K, Cassar G, Borg MK. Mechanical properties of pristine and nanoporous graphene. Mol Simul. 2016;42(18):1502-1511.
[42] Giannopoulos GI, Avntoulla GS. Tensile strength of graphene versus temperature and crack size: Analytical expressions from molecular dynamics simulation data. Proc I Mech Eng Part N: J Nanomaterials, Nanoengineering and Nanosystems. 2017;231(2):67-73.

[43] Brenner DW, Shenderova OA, Harrison JA, Stuart SJ, Ni B, Sinnott SB. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. J Phys: Condens Matt. 2002;14(4):783-802.

[44] Stuart SJ, Tutein AB, Harrison JA. A reactive potential for hydrocarbons with intermolecular interactions. J Chem Phys. 2000;112(14):6472-6486.

[45] Pastewka L, Pou P, Pérez R, Gumbsch P, Moseler M. Describing bond-breaking processes by reactive potentials: Importance of an environment-dependent interaction range. Phys Rev B. 2008;78:161402.

[46] Perriot R, Gu X, Lin Y, Zhakhovsky VV, Oleynik II. Screened environment-dependent reactive empirical bond-order potential for atomistic simulations of carbon materials. Phys Rev B. 2013;88:064101.

[47] Toh CT, Zhang H, Lin J, Mayorov AS, Wang YP, Orofeo CM, et al. Synthesis and properties of free-standing monolayer amorphous carbon. Nature. 2020;577:199-203.

[48] Allen MP, Tildesley DJ. Computer Simulation of Liquids. Oxford University Press; 1987.
[49] Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular Dynamics with Coupling to an External Bath. J Chem Phys. 1984;81(8):3684-3690.

[50] Gamboa A, Vignoles GL, Leyssale JM. On the prediction of graphene's elastic properties with reactive empirical bond order potentials. Carbon. 2015;89:176-87.

[51] Liu F, Ming P, Li J. Ab initio calculation of ideal strength and phonon instability of graphene under tension. Phys Rev B. 2007;76:064120.

[52] Stone AJ, Wales DJ. Theoretical studies of icosahedral C60 and some related species. Chem Phys Lett. 1986;128(5-6):501-503.

[53] Ma J, Alfè D, Michaelides A, Wang E. Stone-Wales defects in graphene and other planar sp2-bonded materials. Phys Rev B. 2009;80:033407.

[54] Leyssale JM, Vignoles GL. A Large-Scale Molecular Dynamics Study of the Divacancy Defect in Graphene. J Phys Chem C. 2014;118:8200-8216.

[55] Buongiorno Nardelli M, Yakobson BI, Bernholc J. Mechanism of strain release in carbon nanotubes. Phys Rev B. 1998;57:R4277-R4280.

[56] Samsonidze GG, Samsonidze GG, Yakobson BI. Energetics of Stone-Wales defects in deformations of monoatomic hexagonal layers. Comput Mater Sci. 2002;23(1):62-72.

[57] Jiang H, Feng XQ, Huang Y, Hwang KC, Wu PD. Defect nucleation in
carbon nanotubes under tension and torsion: Stone-Wales transformation. Comput Methods Appl Mech Eng. 2004;193(30):3419-3429.

[58] Song J, Jiang H, Shi DL, Feng XQ, Huang Y, Yu MF, et al. Stone-Wales transformation: Precursor of fracture in carbon nanotubes. Int J Mech Sci. 2006;48(12):1464-1470.

[59] Los JH, Ghiringhelli LM, Meijer EJ, Fasolino A. Improved long-range reactive bond-order potential for carbon. I. Construction. Phys Rev B. 2005;72:214102.

[60] Zakharchenko KV, Fasolino A, Los JH, Katsnelson MI. Melting of graphene: from two to one dimension. J Phys: Condens Matt. 2011;23(20):202202.

[61] Ganz E, Ganz AB, Yang LM, Dornfeld M. The initial stages of melting of graphene between 4000 K and 6000 K. Phys Chem Chem Phys. 2017;19:3756-3762.