Toughness and strength of nanocrystalline graphene

Pristine monocrystalline graphene is claimed to be the strongest material known with remarkable mechanical and electrical properties. However, graphene made with scalable fabrication techniques is polycrystalline and contains inherent nanoscale line and point defects—grain boundaries and grain-boundary triple junctions—that lead to significant statistical fluctuations in toughness and strength. These fluctuations become particularly pronounced for nanocrystalline graphene where the density of defects is high. Here we use large-scale simulation and continuum modelling to show that the statistical variation in toughness and strength can be understood with ‘weakest-link' statistics. We develop the first statistical theory of toughness in polycrystalline graphene, and elucidate the nanoscale origins of the grain-size dependence of its strength and toughness. Our results should lead to more reliable graphene device design, and provide a framework to interpret experimental results in a broad class of two-dimensional materials.


Supplementary
: Elastic response of nanocrystalline graphene. The stressstrain response of a 2D periodic nanocrystalline samples with grain size 64 Å, loaded at a strain rate of˙ yy = 10 9 s −1 at T = 300K. yy is the applied strain in the y-direction, while σ yy is the measured stress in the y-direction.

Creating nanocrystalline graphene sheets for simulation
We use the following method to generate nanocrystalline graphene sheets with random grain random shapes and orientations. A square sheet of size L and grain size µ has n g = L 2 /µ 2 grains.
First we choose the n g points as 'centers' of the grains at random (distributed uniformly on the sheet of area L 2 ). Then a Voronoi construction with these points is used to generate the granular regions associated with them. Then the orientation of lattice vectors in each grain is chosen at random. The positions of carbon atoms in the grains and at the grain boundaries is assigned with a recently proposed algorithm 1 . This algorithm results in well annealed nanocrystals.

Simulations of nanocrystalline strength
We simulate 24 different combinations of (L, µ,˙ ) (64, 32, 1), (128, 64, 1), (128, 32 The simulations are carried out in the canonical ensemble with constant NVT integration using a Nose/Hoover thermostat with the LAMMPS software 2 . A constant strain rate is imposed in the y direction by using the SLLOD equations of motion 3 . The temperature is set to 300 K.
The interaction between the carbon atoms is modeled by using the AIREBO potential 4-6 , with the modification to the interaction cutoff parameter rc min applied as suggested in Ref. 6. As the applied strain increase, the stress also increases initially. However, eventually fracture is initiated and the sample fails. The peak stress obtained during the loading process is defined as the strength of the nanocrystal. The stress-strain response of a typical polycrystalline sample loaded uniaxially is shown in Supplementary Fig. 3. The thermal component of the stress is subtracted from the net response. Supplementary Fig. 4 shows the loading of a representative nanocrystal in this manner.

Simulations of nanocrystalline toughness
We evaluated the toughness of 500 samples of nanocrystalline graphene each for grain sizes of µ = 16, 32, 64 Å. Nanocrystalline toughness is evaluated by generating a square nanocrystalline graphene sheet of size L = 256 Å with the required grain size and random grain morphology as discussed in the previous sections. A crack tip is introduced at the center of the nanocrystal by applying the deformation field corresponding to a stress intensity factor K I as calculated from linear elastic fracture mechanics, i.e., where r, θ are the polar coordinates of a point with the origin placed at the crack tip, G is the shear modulus, and κ = (3 − η)/(1 + η), and η is the Poisson's ratio. Atoms outside a radius of 100 Å from the crack tip are held fixed at this displacement, while atoms inside this radius are evolved with NVT dynamics. As before, we use the AIREBO potential to model the carboncarbon interaction. The stress intensity factor K I is incremented in steps of 0.1 MPa·m 1/2 , and the system is held at each K I for 1 ps. The critical stress intensity factor, and thus the fracture toughness, is defined as the lowest value of K I for which the crack grows in this manner.
In a recent publication 1 , we grid this space in steps of 0.5 • and generate a grain boundary at each point. We generate GBs with a width of 50 Å. We simulate fracture of the grain boundaries under uniform strain applied perpendicular to the grain boundary. The simulation method is similar to that used for nanocrystalline strength, the difference being that the GBs are periodic only along the GB direction. Thus, to apply a strain loading, a strip of atoms 5 Å wide at the left and right edges of the GB is moved rigidly at the required strain rate. The stress is recorded, an the strength of the GB is defined as the largest stress achieved during the simulation. Since we used a uniform grid, and random grain boundaries sample from a uniform probability density over the (θ 1 , θ 2 ) space, we can readily obtain an approximation for the survival probability, S GB (σ), of a randomly chosen grain boundary. This probability is simply equal to the number of GBs that survived at the stress σ divided by the total number of simulated GBs. A plot of this survival probability can be found in Figure 5 of the main text.

Supplementary Note 2. Theoretical derivations
Detailed derivation of strain-rate and grain-size dependent survival probability of nanocrystalline graphene Consider an isolated defect in the graphene sheet which has a stress dependent energy barrier given by ∆E(σ(t)) associated with it, where σ(t) is the stress at time t. In the case of polycrystalline graphene being considered here, this defect corresponds to individual pentagon-heptagon pairs that make the graphene GBs and TJs as shown in Figures 1b and 2b of the main text. If this barrier can be overcome by thermal fluctuations, then the defect will grow and cause global failure. According to nucleation theory, the probability that the barrier will be overcome in a small time dt is given by ωe −∆E(σ(t))/kT dt, where ω is a prefactor. Thus, the probability that the defect survives during the time dt is given by: Since is the reduced modulus of elasticity for plane stress, Y is the Young's modulus, and η is the Poisson's ration, we can remove the explicit dependence on time from the above equation and get the probability that the defect survives a stress increment of dσ = Y R˙ dt in time increment dt as: Now consider a population of non-interacting defects with a stress dependent density of barrier heights given by f (∆E(σ)), and an area-density given by ρ. There are ρL 2 f (∆E)d∆E defects with barrier height ∆E in a sample of area L 2 . The probability that all of the defects with barrier height ∆E survive the stress increment of dσ, which happens in time increment dt is simply a product of the individual probabilities, and is given by: Thus, the probability that all defects (and hence the graphene sheet) survives the stress increment dσ in time dt is obtained by integrating over the distribution of defects barrier heights, and is given by: Finally, the probability that the sheet survives till stress σ (time t = σ/Y R˙ ) is again obtained by taking a product of the survival probabilities over small increments, which corresponds to integrating the term in the exponential in the above equation, and is given by: To make connection to the theory of extreme value statistics and Weibull distribution, note that the above can be written as: The right hand side of the above equation is of the form F (σ) ρL 2 /˙ , where F (σ) can be identified with the exponential term. An important theorem from extreme value theory 8-10 tells us that for any distribution function F (σ), the following holds under very mild restrictions: where m is a positive real number. Under this assumption (which we verify numerically), we can approximate the survival probability as: where a, b are constants that depend on ρL 2 /˙ .
The scaling suggested in Eq. 3 of the text can be arrived at on the basis of the physical reasoning hinted at in the manuscript; however, here we take a more mathematical approach.
The assumption that Λ(·) is of the Weibull form essentially amounts to assuming that F (σ) has a power law expansion, i.e., then, which gives: Thus, Finally, realizing that the density of defects goes as ρ ≈ c/µ 2 , where µ is the linear grain size, and c is a constant, and normalizing with a reference strain rate˙ 0 gives us: where v = (˙ 0 /cα) 1/m .

Maximum likelihood estimator for survival distribution of strength
Given the survival distribution function for the strength of nanocrystalline graphene S(σ|L, µ,˙ ), the corresponding probability density is: If for a single dataset the data σ i are observed, then the corresponding log-likelihood function is given by i log s(σ i |L, µ,˙ ). If we have q datasets D 1 , D 2 , . . . , D q consisting of n 1 , n 2 , . . . , n q data points, obtained at configurations (L 1 , µ 1 ,˙ 1 ), . . . , (L q , µ q ,˙ q ), then the parameters σ 0 , v, m can be obtained by a joint fit of data obtained by maximizing the following log likelihood function: where the dataset D j consists of observations σ j 1 , σ j 2 , . . . , σ j n j of fracture strengths. The normalization with the number of points in the dataset, n j , is carried out to avoid a bias towards datasets with larger number of observations, because these are typically datasets with smaller values of L that can be simulated at smaller numerical cost.

Maximum likelihood estimator for survival distribution of toughness
The estimation of parameter α for the toughness model (Eq. 3 of the main text) is similar to the discussion in the previous section. However, one key difference should be noted. In the case of strength we were able to get a closed form formula and take the derivative in Supplementary Eq.15 analytically. As the corresponding operation cannot be performed analytically with Eq. 3 of the main text, the derivative has to be taken numerically by using finite differences.