Atomistic aspects of fracture

Any fracture process ultimately involves the rupture of atomic bonds. Processes at the atomic scale therefore critically influence the toughness and overall fracture behavior of materials. Atomistic simulation methods including large-scale molecular dynamics simulations with classical potentials, density functional theory calculations and advanced concurrent multiscale methods have led to new insights e.g. on the role of bond trapping, dynamic effects, crack-microstructure interactions and chemical aspects on the fracture toughness and crack propagation patterns in metals and ceramics. This review focuses on atomistic aspects of fracture in crystalline materials where significant advances have been achieved over the last ten years and provides an outlook on future perspectives for atomistic modelling of fracture.


Introduction
Fracture processes and particularly brittle fracture processes are obvious cases where macroscopic materials properties are almost entirely determined by events at the atomic scale. A propagating crack in a brittle material moves by breaking individual bonds between atoms and can therefore be regarded as a macroscopic probe for the atomic bonding. Nevertheless, textbook analysis of brittle fracture resorts to the thermodynamic equilibrium picture of Griffith (1921), formulating mechanical stability of a crack as a balance between the crack driving force, the energy release rate G, and the surface energy γ s of the two fracture surfaces: G = 2γ s . The driving force on a brittle crack can be obtained from elasticity theory as G = K 2 /E , where E is an appropriate elastic modulus and K is the stress intensity factor characterizing the strength of the stress singularity at the crack tip.
From an atomistic point of view, resistance to crack propagation should be characterised by the forces needed to separate bonds successively at the crack tip. Indeed, the first atomistic studies of fracture (Chang 1970) demonstrated that crack propagation is only possible on certain crystallographic planes within a crystal and that crack tip bonds break after being stretched nearly to their elastic limit (Sinclair and Lawn 1972). Atomistic analysis also showed that the discreteness of the lattice manifests itself in the so-called lattice trapping effect (Thomson et al. 1971). Lattice trapping causes the crack to remain stable until loads K + somewhat larger than the Griffith load K G are reached. The magnitude of the trapping strongly depends on the force law that characterises the atomic interaction (Sinclair 1975) and can explain thermally activated subcritical crack growth (Schoeck 1990;Gumbsch and Cannon 2000). Some success could also be achieved in rationalizing experimentally observed cleavage planes and preferred cleavage directions of the central transition metals based on empirical atomic interaction models (Kohlhoff et al. 1991;Riedle et al. 1996). However, while cleavage can be rationalized, crack tip plasticity in semi-brittle materials can only partially be captured by atomistic models due to the limited length and time scales accessible. These size limitations preclude a full inclusion of a realistic distribution of pre-existing dislocations in the model system. Since the brittle-to-ductile transition (BDT) in semi-brittle materials is largely governed by the relation between loading rate of the crack and the rate of plastic deformation near the crack tip, it obviously relates to dislocation mobility (Roberts et al. 1994;Hartmaier and Gumbsch 2005), which in semi-brittle materials is too sluggish to be reasonably dealt with in atomistic simulations. However, many different processes like dislocation nucleation and dislocation multiplication are involved in the BDT and some of them can be investigated atomistically. For example, atomistic aspects of dislocation-crack tip interaction can be studied for inherently ductile materials like fcc metals, where dislocation motion can be captured in the time scale accessible to MD (Bitzek and Gumbsch 2013).
The study of dislocation nucleation at crack tips (Zhu et al. 2004b;Sen et al. 2010;Yamakov et al. 2014), in particular from crack tip defects like ledges (Thaulow et al. 2011;Gordon et al. 2009), or crackmicrostructure interactions requires very large numbers of atoms to be simulated, and is often modelled with empirical interaction models. The atomic environment near a crack tip, however, deviates strongly from the equilibrium bonding situation. Therefore, many otherwise-reliable empirical potentials fail to correctly describe the fracture behaviour of specific materials, particularly brittle semiconductors and ceramics (Gumbsch 2001;Möller and Bitzek 2014a). On the other hand, such processes and the large systems required are currently well outside the scope of ab initio methods like density functional theory (DFT). Atomistic models which are able to capture such extreme situations have only recently become available (e.g. potentials constructed automatically from QM calculations, either by coarse graining the electronic structure to produce bond order potentials (Mrovec et al. 2007;Margine et al. 2011) or using machine learning approaches (Bartók et al. 2010), and screened potentials which recover the correct bond-breaking behaviour by extending the interaction range as proposed by Pastewka et al. (2012Pastewka et al. ( , 2013, as well as novel multiscale methods that combine DFT-calculations with classical potentials (Bernstein et al. 2009) are increasingly being applied to fracture problems. Such methods allow for material specific simulation studies, i.e., simulations which target the accuracy of DFT calculations and should in principle lead to quantitative results which can be compared to experiments.
A major driving force for such materials-specific atomistic studies stems from recent advances in experiments which have uncovered the existence of crack tip processes whose underlying mechanisms can only be studied by atomistic simulations and of course are specific to the studied material. Examples are the observation of stimulated emission of dislocations from crack tips upon intersection of a lattice dislocation (Michot 2011), the emission of dislocations on glide planes with low resolved shear stresses or not accessible to the crack front (Gally and Argon 2001;George and Michot 1993), and the use of atomic force microscopy both for post mortem analysis, e.g. crack deflection (Kermode et al. 2013) and in situ measurement, e.g. slow crack growth (Célarié et al. 2003;Bonamy et al. 2006;Han et al. 2010;Pallares et al. 2011).
The scope of this review is limited to brittle crystalline materials, which typically cleave to produce atomically sharp surfaces at low fracture energies and in the absence of defects (Lawn 1993). In contrast, fracture of amorphous materials has been reported to proceed by the growth and coalescence of voids, as observed in several AFM studies (Célarié et al. 2003;Bonamy et al. 2006) and in a number of large-scale atomistic simulations with classical interatomic potentials (Ki et al. 2009;Muralidharan et al. 2005). This plastic response is not universally accepted, with recent AFM studies reporting atomically sharp crack tips and no evidence for cavity formation (Guin et al. 2005;López-Cepero et al. 2007). A careful comparison of finite element studies and AFM measurements for silica crack tips showed no evidence for a plastic zone (Fett et al. 2008), and a recent AFM/DIC (digital image correlation) (Hild et al. 2015) study showed that any deviation from linear elastic fracture mechanics (LEFM) is localised to a region within 10 nm of the crack tip (Han et al. 2010). Moreover, the subtle role played by water, whether via deep penetration into the matrix (Lechenault et al. 2011) or by individual crack tip reactions (Michalske and Freiman 1982) is not yet fully understood (Wiederhorn et al. 2013). Fracture of amorphous materials will not be further reviewed here.
Our review focuses on aspects of fracture in crystalline materials where significant progress has been made in the last ten years. We put an emphasis on materials specific simulation rather than generic models for entire materials classes. In the next section advances in the understanding of crack trapping effects will be summarized, thereafter dynamically propagating brittle cracks and the interaction of cracks with other elements of materials microstructure will be reviewed. Finally, chemical aspects of fracture will be reviewed and with a longer outlook we will attempt to point to interesting questions that should be tackled within the next decade.

Bond trapping
Lattice trapping of cracks, introduced above, has been long regarded as the major source of atomistic effects in fracture. Consequences of lattice trapping or bond trapping have been discussed by Gumbsch and Cannon (2000) demonstrating that: (1) lattice trapping leads to an anisotropy with respect to the propagation direction on one cleavage plane; (2) cleavage does not necessarily lead to low-energy fracture surfaces; (3) the ratios of fracture energy to the thermodynamic work of adhesion for perfectly brittle, planar cracks can easily exceed a factor of two especially in non-metals; (4) trapping effects for interfacial cracks can be larger than for crystal cleavage planes.
Cleavage crack propagation in single crystals is known to be an anisotropic process with distinctly preferred propagation directions (George and Michot 1993;Riedle et al. 1996). In tungsten, for example, crack propagation at liquid nitrogen temperature on both the {100} and {110} cleavage planes preferentially occurs in crack propagation directions normal to the 110 crack front (Riedle et al. 1996;Gumbsch 2003). While one may argue that dislocation motion may be associated with cleavage crack propagation in tungsten, cleavage in virtually dislocation free silicon single crystals at room temperature clearly is not associated with dislocation motion. Nevertheless, directional cleavage anisotropy has been observed for both the {111} and the {110} cleavage planes in silicon. The preferred propagation direction is along the 110 directions on both cleavage planes (Michot 1988;George and Michot 1993). This preference has been linked to the low lattice trapping barriers in these particular directions (Pérez and Gumbsch 2000), as calculated by quantum mechanical methods based on DFT. A similar anisotropy was also found in DFT calculations for diamond (Pastewka et al. 2008) and used as a test to validate simpler semi-empirical bond order potentials. On lower symmetry crystals, like potassium bichromate which has an alternating AB structure along {001}, cleavage anisotropy may even be different when propagating a crack forward or backward in one crystallographic direction. Cleavage along a (001) plane in [100] direction leads to upper side being terminated by A while the lower is terminated by B, while cleavage in the opposite direction [100] inverts the surface termination (Plomp et al. 2001).
Most computational fracture studies to date have been limited to narrow, quasi-2D model systems. These models correspond to a straight crack front where all bonds along the front break simultaneously. For dynamic fracture at loads where no lattice-trapping barriers exist, such models can provide very useful insight, but when lattice-trapping effects are significant 2D boundary conditions are clearly unrealistic. An alternative 3D propagation scenario was first suggested by Sinclair (1975) and was later applied in practical calculations for the Si(111)[112] fracture system to calculate the minimum energy pathway (MEP) for crack advance by Zhu et al. (2004a), who showed that the barrier for the formation of a kink pair on a crack front is much lower than the lattice trapping barrier for sequential advance of the entire front. This is because in the latter case the barrier to break a single bond must be multiplied by the width of the system, which is extremely large in macroscopic systems. This kink-based mechanism offers an alternative explanation of the directional anisotropy on the Si(111) plane referred to above, which can be viewed as arising from differences in the energy required to create kink pairs for different crack orientations (Zhu et al. 2004a(Zhu et al. , 2006. A molecular dynamics (MD) based study of fracture in a complex metallic alloy has also observed the formation of kink pairs (Rösch and Trebin 2009). The concept of lattice trapping can be readily expanded to include fracture in environments where the atoms are not arranged on a regular lattice, as is the case at interfaces or in amorphous materials (Schoeck 1990;Gumbsch and Cannon 2000). In such situations the more general term bond trapping is used instead of lattice trapping. Although bond trapping at grain boundaries (GBs) was e.g. assumed to be the underlying reason why polycrystalline oxide ceramics show trans-granular rather than inter-granular fracture at higher crack speeds (Gumbsch and Cannon 2000), only relatively few quantitative calculations of the fracture toughness of GBs in brittle materials can be found in the literature, e.g. Grujicic et al. (1997), Gumbsch (1999), Farkas (2000). Only recently, bond trapping of GB cracks was systematically studied along various large-angle tilt GBs in tungsten (Möller and Bitzek 2014). The authors showed that the fracture toughness for brittle GB cracks depends on the exact location of the crack tip and its propagation direction within the GB-plane (Fig. 1). The local fracture toughness was directly related to the local bonding situation and the symmetry of the structural units of the GB. Under the assumption that the quasistatic GB fracture toughness is determined by the kinetics, e.g. the strongest GB bond, rather than by the thermodynamic equilibrium GB energy, the fracture toughness of GBs was shown to be able to surpass that of single crystals oriented for cleavage on the same crystallographic plane and direction. Whether these findings still hold true for extended, kink-containing GB-cracks still needs to be investigated. Experiments on polycrystalline graphene however show that GBs in 2D can actually show extremely high strengths (Lee et al. 2013).

Dynamical effects-instabilities, defects and perturbations
Undefected crystalline materials sometimes exhibit dynamical effects that can often only be explained by atomistic mechanisms. The material-specific phenomena we focus on here are distinct from well-known high-speed dynamical fracture instabilities, which start above about two-thirds of the Rayleigh speed, and have been observed in a wide range of materials (Ravi-Chandar and Knauss 1984;Fineberg et al. 1991;Marder and Gross 1995). Buehler and Gao (2006) used large-scale atomistic simulations with a generic model material incorporating a transition from linear elastic to nonlinear behaviour to show that nonlinearity (specifically, hyperelasticity) plays a key role in the onset of these high speed instabilities. However, explaining complex, material-specific dynamical phenomena requires a new level of accuracy from atomistic simulations, going beyond the generic features that arise from the discretization of the material lattice (Marder 2004). The strong coupling between lengthscales inherent to fracture means that accurate modelling often requires non-uniform approaches that couple the long-range stress concentration, which can be captured very well by classical interatomic potentials (Singh et al. 2014), with local crack-tip chemistry, which must be modelled with more accurate techniques, typically with quantum mechanical precision (Abraham et al. 1998b;Bernstein and Hess 2003;Kermode et al. 2008). For a detailed review of these "hybrid" multiscale simulations methods for materials problems see Bernstein et al. (2009). Prominent applications of this multiscale approach include the early work of Abraham et al. (1998bAbraham et al. ( , 2000, who demonstrated that a concurrent simulation of dynamic fracture using electronic structure, classical atomistic and continuum finite element methods was possible, and that of Bernstein and Hess (2003), who embedded a tight binding description of the near-tip region within a classical model of a large silicon system to compute lattice trapping barriers.  coupled a ReaxFF (Duin et al. 2001) description of the near-tip region with the Tersoff (1988) interatomic potential to investigate the dependence of crack speed on the available fracture energy in silicon (Buehler et al. 2007).
The 'Learn on the Fly' (LOTF) scheme (Csányi et al. 2004), which uses a similar approach but with crack tip processes typically treated at the DFT level, has been used to study a number of dynamic fracture problems. An early application of the technique showed that nanoscale dynamical instabilities during brittle fracture propagation on {111}-type cleavage planes in single crystal silicon can lead, via a positive feedback mechanism, to experimentally observable surface features (Kermode et al. 2008). The same velocity-dependent instabilities were observed in careful experiments carried out in single crystal silicon specimens, leading to micrometre-scale surface ridges produced by deflection of cracks out of the {111} plane when the crack speed falls below a critical value of about 1,000 ms −1 . A recent analytical solution of the full stress fields arising in this crystallographic orientation and explicitly incorporating the nanoscale crack deflection predicted by the LOTF study also leads to triangular ridges (Chaudhuri 2014). The delicate role played by chemical details of crack tip reconstructions on dynamical crack deflection and path selection has also been investigated using the LOTF approach, explaining the high stability of (111) cleavage above this critical speed (Fernandez-Torre et al. 2010). Kermode et al. (2008) also proposed an atomic-scale explanation for the well-known instability of cracks on the Si(110) surface. Simulations showed that, at low speeds, cracks propagate straight in the [001] direction, while at higher speeds additional energy released by breaking bonds allows the crack to break the slightlystronger bonds required to divert onto inclined (111) planes, which are lower in surface energy than the original (110) cleavage plane.
More recently, the LOTF scheme has been applied to study interactions between propagating cracks and chemical impurities, specifically an isolated substitutional boron dopant in an otherwise-perfect single crystal of silicon. A combined theoretical and experimental investigation showed that chemically induced dynamical crack deflection is possible (Fig. 2), i.e. a single atomic defect can deflect a crack as it travels through a crystal (Kermode et al. 2013). The authors also showed how smooth fracture surfaces can still be obtained when breaking real materials, which inevitably contain defects, by demonstrating that the scattering mechanism is switched off for sufficiently fast-moving cracks.
Our final example of a dynamical fracture phenomenon originating at the atomic scale is the additional dissipative contributions to the fracture energy that arise from the sound waves (phonons) emitted by a  (111) cleavage plane in a silicon crystal containing a single boron defect (orange). Dark blue atoms are treated at the DFT-level, while grey atoms are modelled with the Stillinger-Weber classical interatomic potential (the full model system includes ∼ 10 5 classical atoms, not shown here). Scattering from boron defects leads to velocity-dependent crack deflection (Kermode et al. 2013) propagating crack. Atrash et al. (2011) showed how this can be computed from MD simulations to resolve discrepancies between experimental crack speed measurements in silicon and predictions from continuum elasto-dynamic equations of motion for cracks (Freund 1998).

Crack-microstructure interactions
A quantitative understanding of the relationship between microstructural variables and the fracture toughness of a material is key to improving the performance of both the material and the models describing its mechanical behaviour (Gumbsch 2003;Kumar and Curtin 2007). Cracks interact with second phase particles, other grains, voids or dislocations both through their long-range stress fields and locally with the atomic scale structure of the defects. Whereas the study of elastic interactions between cracks and constituents of the microstructure by continuum methods has become an active field of study (Kumar and Curtin 2007;Belytschko et al. 2009). Only relatively few simulations exist at the atomic level.
Most of the simulations of cracks interacting with a second phase effectively reduce the problem to two dimensions by using periodic boundary conditions (PBCs) along the crack front direction, thus modeling both an infinitely extended crack front and second phase (Rafii-Tabar et al. 2006;Liu and Groh 2014). Two setups are usually studied: either the second phase is situated directly at the crack front or within the path of a propagating crack to study the direct crack-obstacle interaction, or the second phase is located above or below the crack plane, so it interacts with the crack only through changes in the stress field. Although voids are the simplest model for a second phase, only few atomistic studies can be found in the literature, e.g. Abraham et al. (1998a), Liu and Groh (2014). The interaction of a (011)[011] crack in α-Fe with cylindrical voids was recently studied by Liu and Groh (2014), who showed that the location of a void influences the load at which dislocations are nucleated from the crack tip. Machová et al. (2009) studied the interaction of cracks in the same orientation in a thin plate of α-Fe with a rectangular, through-thickness bcc-Cu precipitate. At larger distances, the precipitate did not lead to a shielding of the stress field in front of the crack. However, if dislocations were emitted from the crack tip towards the precipitate it hindered their motion. If the precipitate is located directly at the crack tip, the local stress field is changed and dislocation emission from the crack tip takes place within the precipitate. Zhang et al. (2011) studied the interaction of cracks in α-Fe with fcc γ -Fe as second phase. There, the interphase boundary (IPB) acted as a strong obstacle to dislocations emitted from the crack and dislocation pile-ups lead to the formation of nanocracks at the boundary. None of these simulations however studied cracks on the prefered {100}-type cleavage planes of α-Fe Hribernik (2006), which, in addition to typical problems arising from the accuracy of the embedded atom model (EAM) potential selected (Möller and Bitzek 2014a), further limits the direct applicability of these simulation results to an experimental context.
Other atomistic studies of cracks interacting with second phases relinquished the attempt to be material specific and focused on the study of general features of brittle fracture using fcc metals with cracks oriented such that brittle rather than ductile failure takes place. Petucci et al. (2014) studied the interaction of (001)[100] cracks in Ni with rows of Cu, Pd, Pt, Ag and Au substitutional atoms. The rows were four atoms in diameter, located directly in front of the crack tip and oriented parallel to it, thus again representing a 2D system. In all cases the presence of these four-atom rows significantly increased the critical load to initiate crack propagation compared to pure Ni. This was attributed to the reduction of the peak tensile stress at the crack tip by substitutional atoms. The authors identified the difference in the atomic radius of the substitutional atoms compared to Ni as the central quantity determining the critical load. The stress field around the atomic rows caused by the different size of the atoms was also used to explain the trends in the critical load when the second phase was placed further away from the crack tip. Similar reasoning was also used by Rafii-Tabar et al. (2006) and Musazadeh and Dehghani (2011) to rationalize the interaction of cracks with extended clusters of substitutional atoms. Relating the importance of the difference of atomic radii to other properties of the second phase such as its cohesive energy or elastic constants would require detailed parameter studies. Such studies, which ideally should also include the size and shape of the second phase as well as the properties of its interface to the matrix, are however currently still missing, both for material specific simulations as well as for model materials.
Similarly, the majority of atomistic studies of cracks interacting with second phase particles used quasi-2D geometries. Only very recently, the first simulations have been performed for obstacles that do not extend along the entire crack front. Uhnáková et al. (2014) studied for example the interaction of a 6 nm long (001) [110] crack in α-Fe with a rectangular bcc-Cu precipitate containing 30 atoms. Compared to the situation without precipitate, the precipitate retards the crack propagation. After breaking through the precipitate, crack front waves are visible and the resulting fracture surfaces are rough, in agreement with experiments (Fineberg et al. 2003;Bouchaud 2002). Bitzek and Gumbsch studied the interaction between a propagating 25 nm long crack in Ni and a void with 1 nm radius (Bitzek 2006). Cracks in the studied (110) [001] orientation usually propagate by brittle cleavage, as can also be seen in Fig. 3a,b. However, the interaction with the void locally pins the crack front, leading to a local reorientation of the crack front close to the void. Once the crack front locally attains a 121 orientation, new slip planes become available on which perfectly blunting dislocations can be nucleated, effectively locally inhibiting further crack advance (Fig. 3c). This mechanism leads to a characteristic 'V'-shape of the crack front, which can be found e.g. in experiments on cracks propagating in a temperature gradient in an Si-crystal with the same orientation as used in the simulations (Gally and Argon 2001) (Si has the same slip system geometry as fcc). As this newly discovered mechanism of dislocation emission is only available to propagating cracks, the simulations could explain the differences in the dislocation source configurations observed at crack tips in dislocation-free single-crystal Si (George and Michot 1993;Gally and Argon 2001). In general, very little is know about how the interaction of dynamically propagating cracks with localized obstacles influences the resulting crack surfaces, the crack arrest toughness and the overall fracture behaviour. It is to be expected that atomistic simulations with the possibility of precisely defining the relevant material parameters will contribute significantly to the study of fracture of heterogeneous materials (Ponson 2009).
The interaction of pre-existing cracks with GBs was mostly studied in the context of nanocrystalline metals, see e.g. Farkas (2013) for a recent review. From simulations on pre-cracked 3D nanocrystalline α-Fe samples (Latapie and Farkas 2004), Farkas and coauthors found a combination of intragranular and intergranular fracture independent of temperature and grain size. Crack propagation along the GB was observed if the angle between crack plane and GB was below 45 • , whereas the crack was arrested if the angle was larger than 75 • . In the case of crack arrest, the crack grew under increas-

Fig. 3 Interaction of (110)[001]-oriented cracks in Ni with a void (a-c) and a pre-existing dislocation, (d).
In a-c an initially atomically sharp crack propagates at a constant energy release rate corresponding to 1.1 G c towards a void with radius R = 1 nm. Dislocation emission occurs only when the crack intersects the void and the local crack front orientation is part of a set of glide planes orthogonal to the crack plane (only defect atoms with increased potential energy are shown). d shows the processes which can take place when an initially straight 60 • dislocation cuts the crack front of a stable blunted crack at 0.95 G c (Bitzek 2006;Bitzek and Gumbsch 2013). Color coding is according to common neighbour analysis ing load by nucleation and coalescence of nanovoids on favorably oriented GBs ahead of the crack, similar to the fracture behaviour in fcc metals Farkas et al. (2002).
Only a few authors have studied the interaction of a crack with an inclined GB in a controlled bi-crystal setup, e.g. Miller et al. (1998), Terentyev and Gao (2013). In Al (Miller et al. 1998) as well as α-Fe (Terentyev and Gao 2013), the GBs were shown to act as strong obstacles to crack propagation. This was mostly due to massive plastic deformation that was triggered at the GB by the advancing crack. However, clearly more studies on brittle materials using fully 3D simulation setups and controlled loading conditions as well as a wider selection of GBs need to be performed to evaluate the obstacle strength of GBs with respect to propagating cracks.
On the other hand, numerous studies on intergranular fracture using bi-crystal setups have been published in the last ten years, e.g. Yamakov et al. (2006), Cheng et al. (2010), Paliwal and Cherkaoui (2013) ;Cui and Beom (2014), Möller and Bitzek (2014), Péron-Lührs and Sansoz (2014), including studies on interphase boundary (IPB) fracture (Liu et al. 2013;Yang et al. 2014), see also the article on interfacial fracture in this issue (Banks-Sills 2015). Most of them focus on intergranular fracture in inherently ductile metals like copper (Cheng et al. 2010;Cui and Beom 2014). Due to the different relative orientation of the slip systems in both grains with respect to the crack front, the fracture behaviour depends on the propagation direction of the crack; i.e., whereas brittle intergranular fracture can occur in one direction, the crack can become blunted in the opposite direction, leading to ductile failure only at higher loads. The simulation results can be understood in the framework of the Rice model for dislocation nucleation from crack tips (Rice 1992), which is different from the direction dependent GB fracture toughness in brittle materials caused by bond trapping discussed in Sect. 2 above. However, the additional pos-sibility of dislocation emission from sources within the GB needs to be considered (Cheng et al. 2010;Cui and Beom 2014). The Rice model was also employed to explain the fracture behaviour of a Cu/SiC interface by Yang et al. (2014). Atomistic simulations of GB fracture are also increasingly used to parameterize higher-scale models for interface decohesion, like the cohesive zone model (CZM) (Yamakov et al. 2006;Paliwal and Cherkaoui 2013;Péron-Lührs and Sansoz 2014;Leblond et al. 2015). These models now allow for the direct comparison of finite element (FE) models with fully atomistic simulations of the same nanocrystalline geometry (Coffman et al. 2008;Péron-Lührs and Sansoz 2014). Although the current models are not yet able to capture the entire complexity of GB decohesion in ductile materials, they represent an important step towards the development of a hierarchical multiscale model of fracture. Such models would be of particular interest if they manage to include the effect of segregated atoms on the fracture toughness of GBs and IPBs, see e.g. Farkas et al. (2005), Liu et al. (2013). This would however require a robust link between DFT calculations of traction-seperation data in GBs containing segregated atoms and models of GB fracture toughness to be established, which can be then be used in continuum models, see e.g. Tahir et al. (2013). A common challenge for all meso-or continuum-scale models of GB-fracture is the generalization of the results of very few atomistic simulations, which are usually performed on special Σ-GBs to the entire five-dimensional parameter space characterizing GBs.
Although many authors have studied the nucleation of dislocations from crack tips, the interaction of cracks with pre-existing dislocations has only recently been addressed. Bitzek andGumbsch (2008, 2013) studied the interaction of dislocations on different glide systems with (110)[001]-oriented cracks in a cubic Ni sample with 75 nm side length. For dislocations interacting with static cracks, they were able to identify the following elementary processes which could occur individually or in combination with each other: (a) the cutting of the crack front without further dislocation processes; (b) the cross-slip of dislocations which locally attained screw orientation onto glide planes with higher resolved shear stress, a process which was shown to lead to a spiral source; (c) the partial cross-slip of dislocation segments directly at the crack tip, usually followed by (d) the stimulated emission of other dislocations (cf. Fig. 3d). By performing more than 20 simu-lations of different dislocation crack arrangements and analysing the resolved shear stresses on all dislocations, the authors were able to rationalize which configurations involve which mechanisms (Bitzek and Gumbsch 2013). These simulations helped to explain the stimulated dislocation emission and avalanche-like dislocation multiplication observed in high-temperature experiments on stable cracks in single-crystalline Si at sub-critical load where dislocations were introduced by indentation close to the crack tip (Michot 2011). The situation is different for propagating cracks, where interaction with an orthogonal dislocation leads to the formation of V-sources, similar to the void in Fig. 3a-c (Bitzek and Gumbsch 2013).
The interaction of a propagating (110)[110] crack in Si with a static 90 • partial dislocation was also studied with the LOTF scheme referred to in Sect. 3 above (Makov et al. 2009). Here, the interaction led to a local roughening of the fracture plane, in qualitative agreement with the observations of Shilo et al. (2002). However, the simulations were performed only in a quasi-2D set-up, and the simulations and experiments were performed at temperatures where dislocation mobility in Si is negligible, explaining the differences to Bitzek and Gumbsch (2013).
Recently, Xu and Demkowicz reported simulations in Al that showed that the interaction of cracks with wedge disclinations can lead to healing of the cracks, even under tensile load (Xu and Demkowicz 2013). Whether the rather particular combination of microstructure and boundary conditions necessary for crack healing can also be achieved in experiments remains to be seen.

Chemical aspects of fracture and sub-critical crack growth
While the possibility of chemically activated fracture at subcritical loads has been known for some decades (Lawn 1993;Michalske and Freiman 1982), accurate modelling of crack tip chemistry has only very recently become feasible as a result of advances in both modelling approaches and available supercomputing resources. The effect of environmental species on bond breaking mechanisms is a fundamental issue that is not yet fully understood. It forms part of the more general topic of stress corrosion cracking, a term also used to describe a broad class of complex and interrelated chemical and mechanical processes that can lead to the failure of normally ductile metals and alloys when exposed to an electrochemical environment (Sieradzki and Newman 1985).
Here, we restrict our focus to environmentally assisted subcritical crack growth in brittle materials, typically ceramics such as silica, where chemically aggressive molecules such as oxygen or water react at crack tips (Wiederhorn 1967;Michalske and Freiman 1982;Ciccotti 2009). This is in contrast to other chemical effects leading to embrittlement, such as the segregation of solute atoms at GBs where they reduce the interfacial cohesion, see e.g. Lozovoi et al. (2006) and the discussion in Sect. 4 above, liquid metal embrittlement (Joseph et al. 1999;Luo et al. 2011), or hydrogen embrittlement. In particular for hydrogen, the exact atomistic mechanisms that lead to the reduction in ductility are still under debate (Lu et al. 2001;Wen et al. 2003;Lu and Kaxiras 2005;Song and Wa 2013). We note, however, that there are situations in nominally ductile materials where individual chemical reactions are thought to play a key role in inhibiting dislocation emission leading to embrittlement, for example in aluminium in an oxygen environment (Zamora et al. 2012).
For a long time environmentally induced bond breaking mechanisms were not thought to be important in covalent crystals such as silicon without polar bonds (Cook 2006), despite their vulnerability to fatigue following cyclic loading (Connally and Brown 1992;Muhlstein et al. 2002). Pioneer work simulating crack tip reactions was carried out by Ogata et al. (2001), who developed a hybrid finite-element, MD and density functional theory model and applied it to study the effect of H 2 O on Si (110)-surface crack initiation (Ogata et al. 2004), suggesting that dissociation of water molecules occurs preferentially at crack tips and can lead to Si-Si bond cleavage. However, the relatively small model system and the geometry in which the load was applied complicate a direct comparison with fracture energies measured in experiment. More recently, Colombi Ciacchi et al. (2008) used first principles MD simulations to study how a hydroxylated native oxide layer forms on Si(001) under wet conditions. They determined that water molecules adsorb and dissociate on the oxidised surface, leading to cleavage of Si-O bonds. Interestingly, tensile strain was found to significantly enhance the driving force for dissociative adsorption of water, consistent with a picture where these kinds of reactions lead to environmentally driven subcritical crack growth.
The LOTF scheme discussed briefly in Sect. 3 above was used by Moras et al. (2010) to study how hydrogenfilled "platelet" defects form in silicon crystals following the implantation of hydrogen in the SmartCut TM process used in the semiconductor industry to produce atomically smooth silicon wafers. Dynamical simulations revealed that H 2 molecules form at the internal surfaces of these platelets, diffuse, and dissociate at the ends of the platelets, where they enable stress-corrosive breaking of silicon bonds (Fig. 4a). A similar approach was recently used to demonstrate that cracks in silicon can initiate and propagate when the energy supplied by the mechanical load is insufficient to create new fracture surfaces, providing that oxygen is available (Gleizer et al. 2014). Subcritical crack propagation driven by the dissociative chemisorption of oxygen molecules was predicted by QM-accurate simula- in the presence of oxygen (Gleizer et al. 2014). The O 2 molecule seen in panel (b) spontaneously dissociates, cleaving an Si-Si bond at the crack tip, before chemisorbing on the newly exposed fracture surface as seen in panel (c) tions (Figs. 4b,c) and confirmed by experiments showing subcritical crack initiation in air, but no cracking in an inert, oxygen free environment.
In the technologically relevant case of silica, the nanoscale mechanism for subcritical stress corrosion cracking in an oxygen-or water-rich environment remains a subject of considerable debate (Ciccotti 2009). While stress-driven chemical processes must contribute, it is not known which mechanisms are dominant, and the exact role water plays in crack propagation is not yet fully understood (Wiederhorn et al. 2013). Performing accurate fracture simulations in ionic materials such as silica is complex because of the long range electrostatic interactions. Some progress has been made by considering relatively small model systems, e.g. Zhu et al. (2005) used a semi-empirical molecular orbital approach to show that stress lowers the energy barrier to hydrolysis of Si-O bonds by environmental water molecules in a silica nanorod model system. A parameterisation of the ReaxFF reactive interatomic potential for silica/water (Fogarty et al. 2010) has very recently been applied to model stress-corrosion cracking in silica (Zhang et al. 2014), although comparisons with experimentally relevant loading conditions have not yet been made.
Recently, efficient and accurate polarisable interatomic potentials for oxides have been developed which coarse-grain electronic degrees of freedom, replacing them with dipole moments determined selfconsistently in response to the local electric field (Tangney and Scandolo 2002;Brommer et al. 2010;Kermode et al. 2010). This approach has been applied to model fracture in alumina, where charged fracture surfaces and the piezoelectric response to strain ahead of the crack tip are found to produce complex electric field dynamics (Hocker et al. 2012), which we note could well drive the flow of polar fluids such as water. A full understanding of the stress corrosion problem in silica would need to take such effects into account, as well as the details of the crack tip chemistry and water penetration.

Outlook and future perspectives
The atomistic study of fracture is a highly dynamic field that will continue to profit from advances in both computational power and simulation methods. In addition to continuing work on classical problems like the nucleation and multiplication of dislocations at crack tips or bond-trapping, and on the recent developments identified in this review such as accurate modeling of crack tip chemistry and of the interaction of cracks with increasingly complex microstructures we believe that in the next decade, atomistic simulations will contribute significantly to our understanding of the following scientific problems: 6.1 3D aspects of fracture and effects due to crack front curvature Currently, nearly all simulations are performed on straight crack fronts, often in quasi-2D setups. Unlike LEFM which is an intrinsically two-dimensional theory, this is not a fundamental obstacle but more a problem of computational expense. Crack propagation by kinks (Sinclair and Lawn 1972;Zhu et al. 2004aZhu et al. , 2006 is often neglected in this kind of set-up, which furthermore severely limits the accessible slip systems and the possibility of the crack front to vary its orientation. Crack curvature effects due to the interaction of propagating cracks with localized obstacles or different local rates of crack advance or closure caused by heterogeneities in the stress field have not yet been studied in a systematic way. Crack curvature effects might furthermore become important for the growth of small crack nuclei. Recent studies on penny-shaped cracks (Ersland et al. 2012;Möller and Bitzek 2015) show that dislocations emitted from one part of the crack can interact with those emitted elsewhere as well as with the crack itself in non-trivial ways (Fig. 5). In this context cracks with misoriented (Gordon et al. 2009) or rough crack planes, or defective crack fronts in general are expected to come into the focus of atomistic studies. Although the computational costs are much higher than for quasi-2D setups such simulations could literally open up a new dimension in the study of fracture.

Crack nucleation
Although many scenarios for the nucleation of cracks have been inferred or postulated in the literature (Lawn 1993;Lehmann et al. 2003) and crack nucleation is an inherently atomistic process, only relatively few studies actually try to model the crack nucleation process, mostly at GBs (Dewald and Curtin 2007;Wu et al. 2007;Cheng et al. 2008;Pan and Rupert 2014). In particular at the microscale, many structures can nowadays 24 E. Bitzek et al.

Fig. 5
Plastic events around a penny-shaped crack on the (010)plane in α-Fe. Dislocations emitted from one part of the crack front interact with dislocations or twins nucleated at other parts of the crack as well as with the crack itself (color according to potential energy, only atoms with increased energy are shown) (Möller and Bitzek 2015) be fabricated with nearly no pre-existing crack-like defects or stress-concentrators. Systematic, materialspecific studies of crack nucleation are therefore fundamental for the development of physics-based models of material failure of such structures. Mechanism-based models for crack nucleation are also very important for fatigue (Saxena 2015) in the very high cycle (VHCF) domain, where the fatigue life is almost completely governed by the crack nucleation and initial stages of growth (Mughrabi 2006) and crack nucleation at featureless sites in the interior of the material is often reported (Bach et al. 2014).

Multiscale modeling of fracture
Material-specific quantitative, predictive modelling of fracture requires accurate atomistic models to simulate crack tip chemical processes involving highly stressed bonds as well as inclusion of the microstructural complexity within the plastic zone. This involves the modelling of the motion of many dislocations, their interactions with solute atoms as well as with obstacles like GBs. On the atomic scale, concurrent multiscale models which combine DFT and classical potentials like LOTF (Csányi et al. 2004) or in which continuum models are used to apply complex loading conditions (Kohlhoff et al. 1991;Miller et al. 1998;Möller et al. 2013) have already proven their utility. However, we believe concurrent coupling of atomistic simulations with higher-scale models like dislocation dynamics (DD) or crystal plasticity simulations will only be useful for very specific problems and well defined setups since the time-scale will be entirely set by the atomistic domain. Hierarchical models in which information is passed from one scale to the other seem better adapted to address fracture on all scales. Statistical approaches which use information gleaned from atomistic simulations to coarse-grain e.g. the interaction of a crack front with a field of random obstacles could lead to significant advances in modelling fracture in real, defect containing structural materials. Since failure prediction in real materials is invariably a rare event, proper treatment of statistical considerations is a key feature of a truly predictive multiscale fracture modelling toolkit. At the same time, rigorous consideration of how uncertainties-both due to limitations in model accuracy and to incomplete statistical sampling over possible microstructures-propagate through the multiscale hierarchy is essential.

Materials specific interaction models
Fracture studies to date have mostly relied on generic atomic interaction models (such as EAM potentials) or on material-specific electronic structure studies based on DFT methods. While the latter cannot yet be applied to many of the important fracture problems, the former can not be quantitatively compared with experiments on any specific material. In this context the further development of material specific, chemically accurate but linear scaling atomic interactions using machine learning approaches (Bartók et al. 2010) or bond order potentials (Mrovec et al. 2007Seiser et al. 2013;Cák et al. 2014) are of crucial importance.

Thermally activated processes
In the last years, numerous approaches have been developed to identify minimal potential energy and free energy pathways between states as well as to extend the time scale of atomistic simulations, allowing the direct modelling of thermally activated processes (Voter et al. 2002). Calculation of energy barriers and activation volumes similar to Zhu et al. (2004b) e.g. for dislocation nucleation from defective crack fronts would in combination with transition state theory and statistical approaches be important to model the influence of dislocation nucleation on the BDT. The real strength of accelerated MD methods would, however, most likely lie in the modeling of fracture in amorphous materials, where plastic events with a very broad distribution of activation energies and volumes are involved, e.g. in the formation of voids ahead of a crack tip (where once again there is a clear requirement for realistic 3D model systems).

Importance of dedicated experiments
Although atomistic simulations are expected to contribute to important breakthroughs in the near future, such studies are most convincing if they are closely interlinked with experimental studies dedicated to atomic scale problems. Fundamental advances require the experimental study of model systems like pure single or bicrystals with well-characterized defect content under controlled loading conditions, ideally at a similar scale. Here, the recent advances in in situ micro-and nanomechanical testing in the SEM, AFM or TEM, see e.g. Wondraczek et al. (2006)