Quantitative prediction of the fracture toughness of amorphous carbon from atomic-scale simulations

Fracture is the ultimate source of failure of amorphous carbon (a-C) films, however it is challenging to measure fracture properties of a-C from nano-indentation tests and results of reported experiments are not consistent. Here, we use atomic-scale simulations to make quantitative and mechanistic predictions on fracture of a-C. Systematic large-scale K-field controlled atomic-scale simulations of crack propagation are performed for a-C samples with densities of $\rho=2.5, \, 3.0 \, \text{ and } 3.5~\text{g/cm}^{3}$ created by liquid quenches for a range of quench rates $\dot{T}_q = 10 - 1000~\text{K/ps}$. The simulations show that the crack propagates by nucleation, growth, and coalescence of voids. Distances of $ \approx 1\, \text{nm}$ between nucleated voids result in a brittle-like fracture toughness. We use a crack growth criterion proposed by Drugan, Rice \&Sham to estimate steady-state fracture toughness based on our short crack-length fracture simulations. Fracture toughness values of $2.4-6.0\,\text{MPa}\sqrt{\text{m}}$ for initiation and $3-10\,\text{MPa}\sqrt{\text{m}}$ for the steady-state crack growth are within the experimentally reported range. These findings demonstrate that atomic-scale simulations can provide quantitatively predictive results even for fracture of materials with a ductile crack propagation mechanism.


I. INTRODUCTION
Amorphous carbon (a-C) has many industrial applications, from electrochemical sensors 1 to wear resistant coatings. 2 Mechanical processes, such as plasticity and fracture, play a crucial role in the performance of a-C in these applications. 3 Fracture toughness is particularly important for the reliability of coatings. It is often measured using nanoindentation tests. [4][5][6] While nanoindentation is easy to perform, it is nontrivial to analyze. 7 The values reported for fracture toughness of a-C scatter between K I = 4 − 11 MPa √ m. [4][5][6] Some of the differences can be associated to the difference in density of testing specimen, although this has not been clearly investigated in the literature. The difference between results highlight the fact that it is not straightforward to measure fracture properties using nanoindentation methods.
An alternative to experiments is the calculation of material properties from atomic-scale simulations that are in principle predictive. However, calculating fracture properties [8][9][10] has remained difficult for most materials because the process zone at the crack tip, where the material deforms inelastically, is typically much larger than the accessible simulation cell sizes. Typical sizes range from 1 cm for metals to 100 nm for brittle ceramics while molecular calculations are presently limited to size of few 10 nm's. Molecular calculation of fracture properties for fully developed fracture process zone has therefore been limited to brittle crystals that show little to no plasticity near the crack tip: Silicon [11][12][13][14][15][16][17][18][19][20][21][22][23] and diamond. 24,25 Exceptions are molecular dynamics (MD) calculations on model Lennard-Jones glasses that can yield insights, such as cavitation in front of the crack tip, but are difficult to translate to real-world materials, 26,27 and a few cases of realistic materials such as amorphous Li-Si. 28 Atomic-scale fracture simulations can be divided into multiple categories. Simple equiaxial loading simulations can be used to study bond breaking mechanisms. 29 However, since stress distribution is completely different from the real crack tip, this type of simulation does not reveal the mechanism of fracture at crack tips. The thin strip geometry is often employed for the study of dynamic fracture propagation in MD, due to its simplicity, and the fact that the strain energy release rate is constant and independent of crack length. 20 However, in thin strip loading, the "K-field" near the crack tip has a superimposed uniaxial T-stress which is known to have significant effect on fracture toughness. 30,31 A third solution, which we also employed in the present paper, is to use the solution of linear-elastic fracture mechanics to control the displacement of the outermost atoms of the simulation domain. This approach eliminates the T-stress and can directly provide fracture parameters, given the condition that plasticity is well contained in the simulation zone. This approach goes back to Sinclair 8,32 and was used subsequently to study crystals 13,14,24,25 and amorphous materials. 28 We note that Sinclair 32 used flexible boundary conditions while Refs. 13, 14, 24, 25, and 28 employ rigid boundaries (as we do here). Both techniques yield K c I , but values converge faster with respect to system size for flexible than for rigid boundaries. 33 A systematic study of elasto-plastic behaviour of a-C with different densities was recently performed by Jana et al., 34,35 showing that the stress sustained during steady-state flow of a-C follows a Drucker-Prager law at room temperature. 36 Furthermore, we have also performed uniaxial tests at 0 K to measure the yield stress of a-C. From these results, we estimate that the process zone radius in mode I fracture for a stress intensity factor K I = 10 MPa √ m, which is in the range of experimental results, is less than 22 nm. a-C may therefore be in a sweet spot in which near-tip plasticity occurs but the process zone is small enough such that direct atomistic calculations of early stage fracture are possible. Insights obtained from fracture of a-C may well transfer to other network glasses, such as silica or alumina.
We therefore here apply systematic large-scale, K-field controlled MD simulations to study fracture in a-C as a function of its density. We study initiation and propagation of cracks up to 10's of angstroms of crack length under small-scale yielding conditions, that allow for the measurement of fracture resistance curves and initiation fracture toughness K c I . This simulation needs on the fly calculation of crack-tip position and adjustment of the boundary conditions which is discussed in detail. Furthermore, based on a criteria introduced by Drugan, Rice & Sham, 37 we calculate steady state fracture toughness K ss I under pure mode-I fracture loading. Finally, these calculations enable us to compare our MD results with experimental measurements.

II. SIMULATION METHOD
We use the screened variant of the Tersoff III potential 25,38 for all our simulations. This potential was designed to correctly describe bond-breaking processes 24 which are fundamental to correctly capture plastic deformation and fracture. We perform quasistatic calculations using overdamped MD to simulate fracture in a-C. Quasistatic simulation ensures a clean separation of mechanics from temperature driven relaxation in amorphous systems. On the contrary, room temperature MD simulation of an amorphous material does not capture the real relaxation in the experimental system due to the vast difference between time scales of simulation and experimental settings.

A. Sample preparation
We create our samples by quenching a-C from the melt. The melt equilibrates for 5.0 ps 39 at 6000 K, followed by a linear quench down to 0.1 K at three different rates of 10, 100, 1000 K/ps. Lower quench rates are computationally expensive, thus most of our calculations are done with samples prepared with 1000 K/ps. The slower rates of quench is used to investigate quench rate dependence of the results. Temperature is controlled using the Langevin thermostat during the whole procedure. Finally, we relax the box and atomistic degrees of freedom using energy minimization method. The box is allowed to relax anisotropically with energy change tolerance of 10 −8 eV/atom or root mean square force change tolerance of 10 −8 eVÅ −1 /atom and thus the residual stress is completely removed. Liquid quenches are a standard procedure for the creation of a-C samples, 34,[40][41][42] since direct simulation of a-C thin films growth [43][44][45][46][47] is prohibitively expensive for the sample sizes considered here.

B. Calculation of mechanical properties
Mechanical properties of a-C have been reported in previous works 34,35 at room temperature. However, all of our fracture simulations are performed at zero temperature. Furthermore, it has been shown that the flow stress of a-C de-pends on the internal pressure of the cell which can be built up during quenching and mechanical loading. 35 Specifically Jana et al. 35 observed a Drucker-Prager law with a zero-pressure shear flow stress of 41.2 GPa and an internal friction coefficient of 0.39. Thus the stress is a function of internal pressure and relaxing the internal pressure during the mechanical loading would decrease the measured stresses. In the fracture test we perform here, creation of a crack can lead to pressure relaxation during the test. Thus, to measure the yield stress we apply a uniaxial tension in the z direction and use a fixed periodic boundary in the y direction together with a periodic but relaxed pressure boundary in the x direction. The boundary conditions in the y direction are similar to the fracture test boundary condition that is carried under plane-strain conditions (see Sec.II C) and relaxing the pressure in the x-direction mimics the relaxation that occurs in the fracture test when the crack opens.
Here, we are using zero temperature uniaxial loading to calculate elasto-plastic material parameters of a-C. We use a simulation cell of dimensions l x = 25 Å, l y = 22 Å and l z = 25 Å. This cell is quenched with the method explained in Sec. II A and loaded in z direction by changing the box size in that direction. The box size is relaxed in x direction using a barostat with time constant of 1 ns to remove the stress in that direction. The dynamics is damped with a damping constant of 5 fs.
We will employ the Hutchinson, Rice & Rosengren (HRR) 48,49 solution of elasto-plastic fracture for the interpretation of our results. HRR uses a Ramberg-Osgood 50 plastic material model. We assume that the atomistic samples behave isotropically and J 2 plasticity model is a reasonable approximation of their flow behaviour. Then, the Ramberg-Osgood elasto-plastic material model for a general stress-strain state reads: Here, s i j = σ i j − σ pp δ i j /3 is the deviatoric stress, E is the Young's modulus, ν is Poisson's ratio, n is the hardening power exponent and σ y is the yield strength and we use Einstein summation convention. The parameter α is associated with the yield strain offset. In the Ramberg-Osgood model ασ y /E is the strain after which we consider the material to be plastic. Here we take α = 0.05 which is associated with strain deviation of ≈ 0.2% for ρ = 2.5 g/cm 3 . All other parameters are extracted by fitting the model in Eq. (1) to the simulation results. Relevant extracted parameters are reflected in Table I. The hardening power exponents are first fitted to the data for all samples and calculated to be n 5.5. To have a consistent plasticity model for all samples we fix the hardening power exponent to be exactly n = 5.5 for all densities. Fig. 1 shows the curve to a sample stress-strain test data.

C. Fracture simulations
We study fracture in generalized plane strain conditions. Plane strain creates smaller process (or plastic) zone than   plane stress, thus needs a smaller sample sizes in MD. We use two different simulation cell sizes. The simulation cell is limited by the time required for creating the a-C: Quenching bigger samples is prohibitively time consuming for slower quench rates. The small samples have the dimensions l x = 100 Å, l y = 22 Å and l z = 200 Å, and big samples have the dimensions l x = 500 Å, l y = 22 Å and l z = 300 Å. The cell is periodic in y direction with fixed l y . We put an initial notch in x − z plane with the length of 30 Å along the z direction. Figure 2a shows the geometry and details of the small simulation cell. An example of a propagated crack configuration if shown in Fig 2b. Our molecular calculations follow a technique pioneered by Sinclair 8,32 and subsequently used by many authors. 13,14,24,25,28 We apply displacement boundary conditions on the boundary of the cell compatible with the solution of linear elastic fracture mechanics of isotropic media in mode-I with stress intensity factor K I for a semi-infinite crack field as follows: Here r is the distance from the crack tip and θ is the angle with respect to the plane of the crack. Loading is applied through incremental increases ∆K I of K I , leading to increments of the displacement (∆u x , ∆u z ) on the boundary. The increments of the displacement are applied to the whole sample. Then, the atoms in the boundary region are fixed and the rest are relaxed. The boundary region has a thickness of t = 10 Å in the x − z plane. The atoms in the boundary region are fixed relative to each other and therefore they apply a force to other atoms inside the box. We use increments of ∆K I = 0.05 MPa √ m in each step. The loading step is followed by a relaxation period of 0.5 ps in NVT ensemble with T = 0.1 K and thermostat damping constant of 5 fs. This is followed by an energy minimization step with an allowed root mean square force change of 10 −7 eVÅ −1 /atom for subsequent minimization steps.
The crack starts moving when K I > K c I and the boundary condition should be adapted accordingly. We find the new crack tip position after each step of loading using coordination analysis. Specifically, we detect newly created surface and take the surface atom with the furthest distance along the zdirection as the new crack tip position. Namely, we are taking max i z i where z i are the z-positions of all atoms detected by surface detection algorithm as the new crack tip position.
To distinguish the surface from the bulk, a threshold is calculated based on the mean density of the a-C within an atomcentered augmentation sphere of radius r c = 4 Å. Specifically, given the number of atoms n i within the sphere centered on atom i, the mean density is where m is the mass of a Carbon atom. When atoms are close to a flat surface, a spherical cap with the height h is removed from the respective sphere in the bulk. This reduces the density from the bulk value ρ 0 to ρ surf = (1−ξ 2 /4(3−ξ ))ρ 0 with ξ = h/r c . We define ξ = 0.85 as the threshold parameter for identification of the surface. All atoms with ρ i < ρ surf are identified as surface atoms.
The important parameter in this surface detection algorithm is the cut-off radius r c . It should be taken in the same range as the feature sizes that needs to be detected. A too big r c averages-out the crack-tip features. A small value of r c introduces noise from the quenched disorder of the a-C's structure into the surface detection algorithm. The value of r c = 4 Å and ξ = 0.85 are chosen from trial and error. A similar method for crack tip tracking has been used in previous MD simulation of fracture in amorphous Li-Si. 28 The above method assumes that linear-elastic fracture mechanics is valid some distance away from the crack tip. Analytical and numerical results on elasto-plastic fracture show that this boundary condition remains valid outside plastic zone even for finite deformation elasto-plastic deformation 51 and is valid in nanometer proximity to the crack tip for purely brittle materials. 52 An estimate for the size of the plastic zone under plane strain conditions based on a continuum mechanical analysis is r y = (K I /σ y ) 2 /3π 53,54 in the x-direction for our geometry. Table I shows the anticipated plastic zone size as a function of K I . The maximum possible K I load for each density until which the plastic zone can be accommodated is also indicated in Tab. I for both sizes of samples. Applying the boundary condition of Eq. (2) is realistic for our simulations, i.e. the plastic zone is confined to the region inside the box until the applied load is less than indicated maximum K I in Tab. I.
An additional correction to the continuum solution (which we apply on the boundary), comes from the quenched disorder of the amorphous atomic structure, leading to stress fluctuations throughout the geometry. However, Rice 55 showed that the small deviations from the dominant singular elastic stress, do not have a strong effect on the J-integral or fracture toughness and thus will not alter our simulation results.

A. Simulated R-curves
We used the above-mentioned procedure to measure the crack length a, or rather its change ∆a from the initial condition, as a function of mode-I stress intensity factor K I . This allows us to plot the function K I (∆a), the crack growth resistance curve often named the R-curve. 56 To avoid confusion, one should bear in mind that we are controlling K I and measuring ∆a and the functional notion of R-curve K I (∆a) is used as standard notion in fracture mechanics. A full R-curve manifests three stages of crack growth: (i) The crack initiation regime provides us with K c I . (ii) The crack growth regime is characterized by the slope dK I /d∆a that drops toward zero when the crack enters (iii) the steady state regime where the stress intensity factor saturates to a constant fracture toughness value K ss I . For a semi-infinite crack in a big enough sample loaded according to the elastic K-field, unstable growth only happens if K I = K ss I . In an experimental setting unstable crack growth happens when the rate of change of the energy pumped to the crack is greater than the R-curve slope. This depends on geometry and loading conditions. Furthermore, geometry and loading of the experimental specimen can change the T-stress (the triaxial stress) which significantly changes the steady state fracture toughness. 57 After the crack initiation phase, the crack resistance is a result of energy dissipation in two regions: (i) The fracture process zone, which is a small region close to the crack tip where the fracture processes takes place and, (ii) the plastic zone further away from the crack tip, where plasticity plays the main role and no abrupt change in the micro-structure is visible. The general understanding in the fracture mechanics community is that the amount of energy dissipated in the fracture process zone remains constant and constitutes the work of separation. 56 This is related to the K c I defined above. For ductile materials, the size of the plastic zone increases with the crack advancement. The increase in fracture resistance is mainly due to this contribution. Distinguishing these two effects is the basis for the cohesive zone models of fracture. 56 Due to this separation of effects, it is possible to infer the steady-state fracture toughness K ss I based on the initiation fracture toughness K c I , the initial R-curve derivative dK I /d∆a ∆a=0 and the yield stress and we discuss this point in detail in Sec. IV. Another consequence is that samples with different sizes should lead to the same K c I . Thus K I (∆a) and simulation results before reaching K ss I are useful. We point out that we were unable to reach the steady state in MD simulations due to limitation in computational power which limits the size of the fracture test specimen. Due to these reasons we only continue the fracture simulation up to ∆a ≈ 100 Å and also will be careful in interpreting results of simulations of different sample sizes.
We performed simulations with two sample sizes. The small sample size is chosen in a way to accommodate the plastic zone at the initiation of fracture. The maximum allowable K I load is indicated in Tab. I in the small sample column. The smaller size lets us perform the simulation for very slow quench rates. Figure 3 shows the R-curve of samples with three different densities ρ = 2.5, 3.0 and 3.5 g/cm 3 and three different quench ratesṪ q = 1000, 100, 10 K/ps. Due to the inherent randomness in the material we provide the result of our test for five different samples for each density and quench ratė T = 1000 K/ps that have been melted and quenched independently. The R-curve consists of sections where K I increases but the crack does not advance and ∆a is constant, followed by smaller or bigger jumps of the crack tip position. These jumps are accompanied by the breaking of multiple bonds within the system. The initiation fracture toughness measured in these simulations is smaller than the allowable load for the size of the sample, giving confidence in our methodology. The Rcurve remains linear for higher K I loads which makes it easy to measure initial R-curve slope dK I /d∆a ∆a=0 . Samples with lower quench rates have the same initiation fracture toughness K c I , but tend to have smaller R-curve slope. This effect is however not significant and due to the high computational costs, we do not test slower quench rates for big sample sizes.
The big sample can accommodate plastic zone of larger K I loads which are indicated in Tab. I in the big sample column. Figure 4 shows the R-curve of samples with three different densities ρ = 2.5, 3.0, and 3.5 g/cm 3 . For ρ = 2.5 g/cm 3 , we show all small sample and big sample results in one figure. As can be seen from R-curves for ρ = 2.5 g/cm 3 , initial fracture toughness K c I does not statistically change across all samples. This is in agreement with the separation of effects discussed above. Furthermore, the initial R-curve slope dK I /d∆a ∆a=0 is also the same for small and big samples. Small sample and big sample R-curves deviate after the maximum possible K I load for smaller samples. For visual clarity, we only show the average value of ∆a for each K I for smaller samples and all R-curves for the bigger samples for ρ = 3.0 and 3.5 g/cm 3 . Table II provides K c I and an estimate of dK I /d∆a ∆a=0 . The K c I has been calculated based on the R-curves of the big samples. Due to similarity of the initial R-curve slopes for small and big samples we fit a line to the small sample R-curves with quench rateṪ q = 1000 K/ps and report that as an estimate of dK I /d∆a ∆a=0 . This linear regression is shown in Fig. 3 as a red dashed line. The R-curve slope decreases for the big samples and seems to saturate eventually to a steady state K ss I . Both initial fracture toughness K c I and the initial fracture toughness rate of change dK I /d∆a ∆a=0 increases with the density of sample. Despite variations between different samples with the same density, the increasing trend in K c I and dK I /d∆a ∆a=0 is systematically present.

B. Fracture mechanism
Given the relatively small experimental fracture toughness values for a-C one might expect a brittle fracture mechanism for the material. However, this material can also show a relatively large plastic strain 34 that creates a a condition for the material to undergo a ductile fracture. One important advantage of MD simulation of fracture for a-C is that it enables us to better understand the mechanism of fracture and address questions such as this brittle vs ductile fracture morphology.
Careful examination of mechanism in MD simulation of a-C shows a consistent mechanism for all densities. A nanosized void is formed in front of the crack tip and then crack grows by coalescence mechanism. Figure 5 shows the void formation for three densities in the course of crack growth. It is clearly visible that ahead of the crack tip there is a sharply localized zone in which the mean coordination number drops as an indication of void creation. The small size and distance of voids from the crack tip is responsible for low fracture toughness values and brittle like fracture (see Sec. IV). The size and the distance of voids from crack tip only slightly changes with density, but generally increases with density.

A. R-curve trends and qualitative analysis
Performing careful simulation of mode I fracture loading constitutes a reliable simulation of the fracture process zone. When we specify an elastic fracture zone with our choice of boundary conditions, a plastic zone emerges inside the sample and a process zone forms near the crack tip. The experimental value for fracture toughness of a-C is rather low, which is consistent with our simulations. Thus one expects brittle fracture to take place. On the contrary, void creation, growth and coalescence mechanism which we report for this material are the typical features of ductile fracture. According to the Gurson model for ductile fracture, the initiation fracture energy (akin to initiation fracture toughness K c I ) is proportional to void spacing l * and yield stress, i.e. J c I ∝ σ y l * . 58 In our material there are no preexisting voids, but voids are created in the fracture process zone at nanometer-scale distances from the crack tip. This void creation process suggests a nanoscale value for l * . Thus, the Gurson model predicts very small values for K c I which resembles brittle fracture despite the ductile fracture mechanisms.
Qualitative trends in our R-curves are consistent with the behavior expected from material properties. From the Gurson model, we expect bigger values for the initiation fracture toughness K c I for higher densities: The yield stress and void spacing l * increase with density and thus according to the Gurson model, we expect an increase in the initiation fracture toughness with increasing density. Furthermore, the increase in dK I /d∆a ∆a=0 for higher densities is consistent with the increase in yield stress and thus higher dissipated energy in plastic zone, and also higher energy needed for tearing of material in the fracture process zone. The material with higher  yield stress and similar micro-mechanism of fracture dissipates more energy in the plastic zone during advancement of the crack.

B. Quantitative predictions and comparison to experiment
All of our simulations are carried out quasi-statically, essentially corresponding to the 0 K limit. There are two competing effects in the mechanical behaviour of amorphous a-C, namely the fast quenching rate which decreases the yield stress and the zero temperature effect which tends to increase the yield stress of a-C. These competing effects may imply that our simulations on rapidly quenched a-C at zero temperature are more relevant than one might envision. Nevertheless, we will be careful in correlating fracture and yield stress to extract trends that are related to experimental results. Furthermore, room temperature simulations would not contribute to a better understanding of the phenomena due to the fact that the competing effects of temperature relaxation and strain rate are occurring at completely different scales in MD simulations and experiments.
There are three reasons that make it difficult to make a direct comparison between MD calculated fracture toughness and experimentally measured values, namely: (i) different geometry of experimental samples which creates complex loading conditions (ii) T-stresses in experimental setting and (iii) different yield stress in MD samples. In particular, a positive T-stress significantly reduces the fracture energy. 59 While very high values of positive residual stress has been reported in deposited a-C thin films, 60 this is not taken into account for fracture toughness estimates from nano-indentation results. In general, we expect a deviation from experiment due to Tstresses and geometrical factors. To be consistent, and rule out the geometrical factors, the two values of K c I and K ss I provide the lower and upper bounds for experimentally measured fracture toughness independent of geometry with zero T-stress. The values for K c I shown in table II for different densities are all in the lower range of experimental values. Given the fact that the effect of density on the fracture of a-C has not been investigated experimentally, the simulated values provide us with some insight about the effect of density on initiation fracture toughness.
In order to find the steady state fracture toughness, we use a criterion introduced by Drugan, Rice & Sham 37 for crack propagation which is similar to crack growth criterion based on the attainment of a critical accumulated plastic strain. 37 Using this criterion and small scale yielding assumption the crack growth follows a nonlinear differential equation: 37,61 Here T is defined as tearing modulus and T 0 is the initial tearing modulus T 0 = T ∆a=0 . Values of α and β depend on the material properties, but α/β ≈ 0.1 is a good approximation for ν ≈ 0.25. 61 We can find the steady state fracture toughness, using the limit of Eq. (3) when ∆a approaches infinity,  Fig. 4. In Fig. 4, we also show the upper and lower bounds of the model as shaded area. The shaded area shows the 95% confidence interval using two standard deviation ±2σ difference in K c I and dK I /d∆a ∆a=0 (reflected in Tab. II).
The values estimated in Table III for K ss I capture the range of reported experimental values in the literature very well, i.e. K I = 4 − 11 MPa √ m. [4][5][6] This agreement suggests that the fracture of a-C is very localized and the MD cell sizes used here can capture all the important effects during crack growth. Table III shows that the extrapolated values for K ss I are lower than the calculated R-curves in Fig. 3. This may be related to the model assumptions of Eq. (4) or ambiguities in the measurement of the yield strength for our atomic-scale a-C models. For example, to the best of our knowledge, there is no direct experimental measurement of the yield strength of a-C in the literature. On the other hand, calculation of yield stress in MD is not straight forward (explained in Sec. II B). If we take the fixed volume condition for yield stress measurement (similar to Jana et al. 35 ) the predicted values of yield stress increases over our estimates and thus K ss I would decrease according to Eq. (4). In contrast, performing a test by completely relaxing the pressure during the test would lower the yield stress and thus increase the estimated K ss I according to Eq. (4). Here we choose a trade-off between these two extreme cases by fixing the l y and relaxing l x which we believe is the best predictor of the yield stress in the fracture test.

V. CONCLUSION
In this paper we performed large scale MD mode-I fracture simulation of a-C with K-field controlled boundary conditions and a well-confined plastic zone. Since experimen-tal measurements of fracture toughness in a-C as a function of density remain challenging, our simulation can provide valuable insight about the fracture process. In this work we simulated fracture of a-C with three different densities ρ = 2.5, 3.0 and 3.5 g/cm 3 . The fracture mechanism for all samples is void creation, growth and coalescence. However, due to the small distance between voids, the fracture toughness remains in the usual range for brittle fracture. Both calculated K c I and extrapolated K ss I are in the range of experimentally measured values which is an indication that the fracture process depends on localized effects near the crack tip that are captured in our atomic-scale simulations.
We want to emphasis that the quantitative predictive power of the atomic-scale simulation results of the fracture of a-C is very encouraging and the methodology and insights are applicable for other glassy materials with similar properties. Other than that, atomic-scale simulations allow controlled testing of various assumptions underlying fracture mechanics, e.g. one could test mixed mode fracture scenarios, or test different fracture path selection criteria such as maximum energy release rate or principle of local symmetry. We also propose more accurate experimental methods, such as micro beam bending experiment that have been performed on hydrogenated a-C, 60 to better understand fracture of a-C.