Ductile-brittle transition at edge cracks (001)[110] in BCC iron under different loading rates in mode I: a 3D atomistic study

3D atomistic simulations via molecular dynamics (MD) at temperature of 0 K and 295 K (22 °C) with a high quasi-static loading rate dP/dt of 2.92 kN s−1 show that cleavage fracture is supported by surface emission of oblique dislocations 111¯011 and by their subsequent cross slip to {112} planes, which increases separation of the (001) cleavage planes inside the crystal. Under the slower loading rate by a factor 5, the crack growth is hindered by twin generation on oblique planes {112} and the fracture is ductile. The MD results explain the contribution of the crack itself to the ductile-brittle transition observed in our fracture experiments on Fe-3wt%Si single crystals of the same orientation and geometry, loaded at the same rates dP/dt as in MD. The loading rates are equivalent to the cross head speed of 5 mm min−1 and 1 mm min−1 used in the experiment. The MD results also agree with the stress analysis performed by the anisotropic LFM and comply with experimental observations.


Introduction
The ductile-brittle transition at cracks in bcc iron-based materials caused by temperature or by loading (strain) rate is often studied via experiments and also theoretically, because the topic is extremely important for engineering applications [1]. A broad overview concerning the influence of loading or strain rate on material properties can be found e.g. in [2], including the effect of loading rate on ductile-brittle transition and fracture toughness.
This atomistic study was motivated by our fracture experiments [3,4] on edge notched single crystals Fe-3wt%Si (see figure 1) loaded in mode I with various cross head speeds from 0.5 mm min −1 up to 5 mm min −1 , where the ductile-brittle transition was monitored at room temperature with increasing loading rate for crack orientation (001) [110] (crack plane/crack front)-see figure 2. The crack plane (001) is the most frequent cleavage plane observed in bcc iron [1]. In particular, in the first experiment [3] with crystals of length L≈100 mm, an unexpected sudden brittle fracture at Griffith level was detected in the case of the longer edge crack (001) [110] with the relative ratio a/W=0.42 (figure 1) under the highest cross head speed of 5 mm min −1 . The presented new free 3D atomistic simulations using molecular dynamics (MD) technique are focused on the experiment in [3]. The simulations utilize the self-similar continuum concept from linear fracture mechanics (LFM), valid for both experimental (macro) specimens and atomistic samples (of sufficiently large dimensions), since in LFM the stress concentration depends just on relative sample dimensions such as a/W, L/W and B/W in figure 1. The previous plane strain MD simulations presented in [3] were of 2D character, where the atomic motion in x 3 =[110] direction was not allowed and the dislocation emission on the slip systems oblique to the crack front was not enabled. 2D simulations in [3] resulted in cleavage fracture, accompanied just with some attempts to emit dislocations on the inclined slip systems, containing the crack front. Here we present new 3D simulations unrestricted by plane strain conditions along the crack front, i.e. with free surfaces like in experiments. Our question is how the crack itself can contribute either to ductile or to brittle behavior in dependence on the loading rate in 3D. In this paper we bring new information, presented in section 3 and briefly mentioned in the Abstract.   [3,4] with edge cracks (001) [110] loaded in mode I at room temperature: (a) Kcr versus VE (mm min −1 ); (b) Kcr versus strain rate (min) −1 . The black points come from [3], where longer crystals of the length L∼100 mm were used. The pink points come from [4], where smaller crystals of the length L∼50mm were used. The transition is indicated by the dashed line.
On the other hand, a direct quantitative comparison between MD results and fracture experiments is not possible due to different strain rates, pre-existing dislocations, impurities and other defects, present in the experiments. A review paper on the capabilities and limitations of MD simulations can be found in [12]. In this paper we treat a difficult crack orientation (001) [110] that is sensitive to strain rate in fracture experiments [3,4]-see figure 2.
As for strain rates in the classical fracture tests with quasistatic continuous loading, the monitored P -t dependencies (resulting force P versus time t) show that brittle fracture occurs within a time interval of several seconds-see e.g. figure 7 in [3] or figure 5 in [4] in the case of iron crystals. Since MD simulations need to use a very fine time step h of the order 10 -14 second to secure numerical stability of Newtonian equations of motion in the atomic lattice (with interatomic distances of the order of 10 -10 m), modelling of real fracture experiments via MD is impossible because of high time consumption (10 14 time steps h are needed). On the other hand, the continuum theory of elasto-dynamic crack problems [13] states that a fast or slow loading is a relative term, depending on the ratio of the critical rise time Δt and the half period T λ/2 of the basic sample vibrations. T λ/2 in mode I depends on the sample length L (λ/2=L). For Δt ? T λ/2 , the influence of inertial forces in Newtonian equation of motion is relatively small and the loading can be considered as a quasi-static loading and vice versa, for Δt = T λ/2 , it represents an impact loading. The validity of this idea can be verified by comparing MD results with the continuum LFM solutions in two limit cases: i/ for impact loading (solved before in [14]), and ii/ for the static or quasi-static case (solved before in [6] and also presented here).
To answer our question concerning the effect of loading rate on the ductile-brittle behavior of edge cracks (001)[110], we utilize in MD the elasto-dynamic concept described above and introduce a time similarity with the fracture experiment from [3] via the relation = dP dt dP dt ex MD where index 'ex' refers to the experiment and MD to the 3D atomistic simulations. In MD, two different quasi-static loading rates 1 PH and 5 PH are treated, corresponding to cross head speed of 1 mm min −1 and 5 mm min −1 in the experiment. MD simulations were performed both at temperature of∼0 K to enable a comparison with continuum solutions and at room temperature of∼295 K (22°C) to enable a comparison with experiment. Further details of MD simulations are described in section 2. The results in section 3 are discussed not only in relation to the mentioned experiment [3] and to its further post-fracture experimental analysis via transmission electron microscopy (TEM) and scanning electron microscopy (SEM) in [15] but also in relation to the recent experiments on Fe-Si single crystals presented in [16]. [¯] (the width W in figure 1) and 30 planes (110) along the crack front, parallel to the direction x 3 =[110] and the sample thickness B in figure 1. The total number of atoms in the atomistic sample is N=2 025 000. The pre-existing edge-through crack of length a was created by cutting the interatomic bonds in the [001] direction along the crack faces (i.e. between layers 450 and 451). During the simulations, the interactions across the initial crack planes (001) were not allowed. The relative dimensions of the atomistic sample a/W, B/W and L/W correspond to the relative dimensions of the specimens treated in fracture experiments in [3], so that the boundary corrections factors F I (a/W) and f I (a/W) staying in LFM at the stress intensity factor K I [17] and the T-stress [18,19] are the same in MD and experiment. Note that F I and f I in LFM do not depend on sample length if

MD simulations
which is obeyed in this study. As follows from continuum predictions in [11], a reliable description of thermal activation at dislocation emission from the crack tip in MD requires a sample thickness / is the interplanar distance in the [110] direction along the crack front. The condition above is obeyed and the used thickness B is sufficient for a reliable description of the thermal activation, which is presented in [8] for temperature of 300 K and 600 K.
When a 3D crystal with free surfaces is generated, unresolved interatomic forces occur at the free surfaces due to missing interatomic interactions. Before loading or heating, surface relaxation in the atomistic 3D sample is performed to bring it to equilibrium and to avoid the influence of surface relaxation on the crack tip processes in subsequent simulations. This was done using the viscous damping in Newtonian equations of motion. The subsequent simulations are free, i.e. no damping is used in Newtonian equations of motion for the individual atoms.
Newtonian equations of motion for the individual atoms are solved in all 3 directions x 1 , x 2 x 3 , using a stable time integration step h=1×10 -14 s, tested in [8].
For the description of the interatomic interactions an N-body potential from table 2 in [20] is used with the lattice parameter ao=2.8665 Å and the basic elastic constants C 11 =243 GPa, C 12 =145 GPa and C 44 =116 GPa for bcc iron. The relevant surface formation energy with this potential is γ 001 =1.812 J m −2 . The corresponding values of Griffith stress intensity factors calculated from the relation g = CK 2 G 001 2 are K G =0.9115 MPa m 1/2 for plane strain and K G =0.835 MPa m 1/2 for plane stress conditions in the atomistic samples. The anisotropic elastic coefficient C depends on crystal orientation and stress state: C=5.198×10 -12 m 2 /N for plane stress and C=4.362×10 -12 m 2 /N for plane strain in our case (for more details, see e.g. [6]).
The relaxed atomistic 3D sample is loaded either at temperature of∼0 K (i.e. without previous heating) or at room temperature of∼295 K (22°C), i.e. with previous heating to this temperature before loading.
As mentioned in Introduction, the different loading rates in MD are introduced via the relation = dP dt dP dt .

ex MD
This is calculated with Young modulus E =175.5 GPa, B∼2 mm, W∼10 mm, L∼100 mm and for the cross head speed VE=dV/dt either 5 mm min −1 or 1 mm min −1 as in the experiment [3]-see figure 2. This elastic approach corresponds to either dP/dt=2.92 kN s −1 (denoted as 5 PH for 5 mm min −1 ) or to dP/dt=0.584 kN s −1 (denoted as 1 PH for 1 mm min −1 ). These values are considered to determine the loading rate (dσ/dt) in MD via relations ( ) ( ) / to each atom li in the 6 upper and lower BW surface layers (001) in figure 1. This distribution of external forces enables us to avoid damage at the loaded BW borders.
In thermal MD simulations, the relaxed atomic configuration is heated to an initial temperature T i by prescribing random atomic velocities to individual atoms according to Maxwell-Boltzmann distribution. The initial kinetic energy in the ensemble corresponds to Nk T 3 2 , where k B is the Boltzmann constant and N stands for the total number of atoms given above. During heating period, the system is not loaded, i.e. Newtonian equations are solved with zero external forces. During the time integration of Newtonian equations, the total kinetic energy Ekin=Ekin(t) and the potential energy Epot=Epot (t, 0)=Epot(t) -Epot(0) in the system fluctuate [8]. With increasing time, Ekin gradually decreases and Epot increases since the atoms deviate from their equilibrium positions and the atomic lattice is subject to thermal expansion. After some relaxation time, a steady state is reached in the system, corresponding to desired temperature T. Here T is given by Boltzmann relation for three degree of freedom á ñ = Ekin t Nk T 3 2 , where á ñ Ekin t ( ) is the time averaged value of the total kinetic energy in the system in the steady state (see e.g. Figure 1(b) in [8]) with á ñ Ekin t ( ) =1.237×10 -14 J and T=295.1 K at time t=2000 h. In statistical mechanics (thermodynamics), the heated atomistic sample in the steady state corresponds to canonical ensemble (N,V,T) with a constant number of particles N, constant temperature T and constant volume a )where α is a coefficient of thermal expansion, which for the used potential agree well with experimental data for bcc iron [21]. This character of the ensemble remains practically unchanged after loading, as presented and discussed in [8].
After heating th edge crack in MD undergoes thermal vibrations. We decided to start with external loading of the 3D atomistic sample at time step 2000 where temperature corresponds to temperature of∼295 K (22°C) and the monitored maximum crack opening displacement (COD) fluctuates around zero values. (The tests with a smaller time step of 0.5×10 -14 s did not improve the results.) Beside COD, also energy balance in the system, i.e. Epot, Ekin and the work Wext done by external forces are monitored each time step, further the crack length in the middle of the crystal (at B/2, layer 16) and the total number of interatomic interactions in the system LINT.

Results and discussions
Except for the core zone where singularity occurs at crack tips or dislocation cores in continuum models, the continuum solutions should be valid at further distances from the core on a macro-scale and also on the atomic scale under equivalent boundary conditions, namely in the case of self-similar linear fracture mechanics (LFM). The comparison of LFM and MD solutions at further distances from the core zone allows us to verify the atomistic and continuum models mutually.
In section 3.1, we present the results of MD simulations at temperature of∼0 K focused mainly on the comparison with continuum predictions from LFM, that do not take into account the random thermal atomic motion. The results at∼0 K from section 3.1 are also needed for comparison later in the section 3.2 to assess the effect of thermal activation at room temperature. The comparison of the stress field by LFM and MD was performed at applied stress level of 1 GPa. The applied stress is within the linear (elastic) region of the theoretical stress-strain dependence for the potential used, presented in figure 1 in [14]. The faster linear loading 5 PH (2.92 kN s −1 ) showed no change in the total number of interactions LINT monitored via MD code. Stress calculations on the atomistic level were performed using the inter-planar concept presented in [22]. We stopped the linear loading at time step 6332 h at the level of σ A =1 GPa and using the critical viscous damping in Newtonian equations of motion we obtained the static atomic configuration in the 3D atomistic sample loaded in mode I. Comparison of MD static stress components S x , S y and T-stress on the x-axis (angle θ=0) with the continuum predictions by anisotropic LFM [6] for the (001) and the function Re(μ 1 μ 2 ) of the complex variables μ 1, μ 2 depends on plane strain or plane stress state in the atomistic sample. Figure 3(c) shows the instantaneous stress component S33 from MD on the x-axis at the applied stress of 1 GPa.
Here it is obvious that plane strain state with σ z >0 prevails only in the nearest vicinity of the crack tip where the peak occurs, while in the larger part of the sample plane stress state prevails with σ z ∼0. For that reason we consider m m = -Re 1.2868 1 2 ( ) valid for plane stress [35]. The boundary correction factor is = [18,19]. The LFM solution is denoted in figures 3(a) and (b) by lines, MD results via circles (o) and r means the distance from the crack tip on the x-axis, oriented in the 110 [¯] direction. Regarding MD, some initial stress arises at the crack front in the atomistic sample after the surface relaxation. The maximum values at the crack tip of the longer crack are = S10 2.266GPa, = -S20 0.863GPa, = - ) Relation for T MD follows from equation (1). The agreement between the static atomic stress S x , S y and the LFM components σ x , σ y is very good in the interval = á ñ r do 1, 10 . / Here r do / is a dimensionless distance where = do ao 2.As for T-stress, in the nearest vicinity of the crack tip, an anisotropic solution s m m = +S S T Re A x y 1 2 ( ( ) ) / for the biaxial loading from [6] rather prevails, while the isotropic static solution prevails with increasing distance r/do from the crack tip toward the right free surface BL. The agreement mentioned above is not trivial; it confirms the self-similar LFM character (no information on boundary correction factors exists in MD simulations) and justifies the concept used in this study. Under the linear loading of 5 PH (2.92 kN s −1 ) up to a high stress level the crack initiation in MD was monitored already in the vicinity of the lower Griffith stress 1.84 GPa, corresponding to the Griffith stress intensity K G =0.835MPa m 1/2 for plane stress, i.e. at s =1.84 GPa which follows from LFM for cleavage fracture. To investigate the microscopic processes at crack initiation and to avoid discussions of dynamic effects, such as flight time corrections of the loading waves, etc, we use a quasi-static ramp loading presented in figure 13 for temperature of 295 K. In linear region, the loading rate corresponds to 5 PH (2.92 kN s −1 ), which means a stress rate 1.84 GPa/11650 h for our atomistic 3D specimen. At temperature of 0 K, the loading starts at t=0 h. For t11650 h the applied stress in MD is kept at a constant level SA=1.84 GPa. Brittle fracture in this run is illustrated in figure 4 by overall shape of the fractured atomistic sample at the end of the simulation. The dependencies RK-time step and RK-COD in figures 5(a) and (b) show that the sample was fractured at a time step∼17650 h (when RK∼0) and when COD corresponds to∼40 Å. The peak value RK max =0.336×10 -6 N corresponds to a global critical stress in the atomistic sample s = RK BW MD max ( ) / =1.82 GPa, which agrees very well with the independent LFM prediction σ G =1.84 GPa for cleavage fracture. This means that the peak value RK max can be related to fracture toughness of the atomistic sample via the relation .825 MPa m 1/2 and G MD =3.54 J m −2 using the elastic stiffness coefficient C=5.198×10 -12 m 2 /N for plane stress. These data are close to Griffith values K G =0.835MPa m 1/2 and 2γ 001 =3.612 J m −2 expected for cleavage fracture. Nevertheless, graphical treatment of MD results shows that this is not a pure cleavage fracture.
The crack was initiated in the middle layer 16 (figure 6(b)) already at time step t i =10978 h, while no crack advance is monitored in the first surface layer 1 in the [110] direction due to dislocation emission on á ñ 111 011 { } oblique slip systems, shown in figure 6(a). The situation in the last surface layer 30 is similar to layer 1. A scheme for this emission is in the right part of figure 6(a). The emission can be thought like a block like shear (BLS) process (see [23,24]), where the upper plane (011) containing the blue crack tip atom A is rigid and the lower plane (011), containing the red atom B, is the moving plane in the direction of the Burgers vector b oriented in the 111 [¯] direction. 3D visualization of the dislocation emission from the edge crack is presented in figure 6(c) using the local numbers of interactions KNT(li) of the individual atoms li from MD with the used short ranged potential [20]. Dislocation cores in figure 6(c) are shown via the red atoms with KNT=16, in agreement with BLS simulation of dislocation generation in perfect crystals. The atoms lying on free {110} surfaces with KNT=10 are blue and the atoms lying on the free crack surfaces (001) with KNT=9 are yellow. Note that the KNT numbers are valid just for the short ranged potentials, where the coordination number at one atom in the perfect bcc lattice is KNT=14, as in our case. The coordination numbers KNT at the atoms lying on free surfaces are lower due to missing interatomic interactions.
As it follows from the geometry of the bcc lattice, the nonzero relative displacements at dislocation emission in figure 6(a) are the components D = - . The relative displacement of the lower part of the crystal with respect to the upper part is { } In agreement with Rice model [9], we consider the emission to be complete when the relative shear displacement U10 in the 111 which is shown in figure 7 by the horizontal line at the level 1. Recognition of slip patterns in 3D can be performed in small perfect bcc crystal via independent block like shear (BLS) simulation mentioned above. The BLS simulations are described in detail e.g. in [24] and presented here in figure 8. Figure 8 [¯] it creates the pattern shown in the middle in figure 8(a). This is the same pattern as the lower oblique slip trace in figure 6(a) after dislocation emission on the oblique plane (011). The upper slip patterns in figure 6(a) come from the mirror oblique plane 011 .
(¯) The oblique dislocation is initially emitted from the surface crack tip as a short curved dislocation that extends gradually under the shear stress on the (011) plane. When such dislocation reaches the straight screw character, parallel with the Burgers vector b= 111 [¯] shown in figure 8 (this vector lies both in the (011) and (112) slip planes), the screw dislocation may change the slip plane to (112) via the cross slip. BLS patterns arising in our observation plane (110) after the slip 111 112 [¯]( ) are shown in figure 8(b). This slip creates a horizontal trace given by the horizontal intersection line of planes [¯] direction of the crack extension. Figure 8(b) illustrates that the slip leads to separation of planes (001), which can help with cleavage along these planes. This is illustrated in figure 9 for the cracked MD sample at a later time 13000 h after the crack initiation and multiple dislocation emissions at the free WL surfaces (110). The multiple surface emissions to oblique á ñ 111 011 { }slip systems and the subsequent cross slip of the emitted dislocations to oblique slip planes {112} in figures 9(a) and (c) causes at the interior cleavage planes (001) the mentioned separation, which assist the crack extension in the interior of the crystalsee the detail in the middle layer 16 along the crack front in figure 9(b). Again, the scheme in figure 9(b) is valid for the lower slip patterns. The upper horizontal slip patterns come from the mirror oblique plane 112 .
(¯) The process described above is a new (unpublished) mechanism for a brittle extension of the edge crack (001) [110].
The process is possible also in anisotropic LFM at Griffith level K G , where K G is considered a material parameter that does not depend on sample dimensions. In the oblique slip system 111 011 [¯]( ) slip stress τ b can be obtained by transformation of the stress tensor from the original coordinate system in figure 1: x 1      [¯]( ) ( ) / /  { } (b) middle layer 16cross slip contributes to separation of the (001) planes, which facilitates fracture; (c) situation in the last surface layer 30 is opposite to layer 1.

In slip system 111 112
[¯]( ) the Schmid factor is∼0.47, thus it is larger than in the original slip system 111 011 [¯]( ) where Schmid factor is∼0.41. This is a qualitative explanation why the cross slip is realized in MD. Figure 6(a) shows that the cross slip begins at the distance of about r=8b from the crack tip under the critical level of loading K I =K G =0.835MPa m 1/2 , which gives by LFM relations (4) and (5): τ b =3.51 GPa. This value is larger than Peierls stress τ P ∼2 GPa needed for motion of screw dislocation in a perfect bcc Fe crystal at temperature of∼0 K, as presented in a recent comparative study [25]. On the other hand, the LFM relation (5) shows, that for the distance r=b/2 from the crack tip and the same loading K G =0.835MPa m 1/2 , the stress concentration in the 111 112 [¯]( ) slip system is lower (τ b ∼14 GPa) than in the slip system 111 011 [¯]( ) treated above. Moreover, the value of∼14 GPa is not sufficient for a direct generation of dislocation in the slip system 〈111〉{112} where the stress barrier is τ disl =16.3 GPa [6] with the used potential. In other words, the anisotropy of the stress field described by LFM relations (3), (4) and (5) explains why the surface dislocations are initially generated in the oblique slip systems á ñ 111 011 { }and later perform the cross slip to oblique system á ñ 111 112 .
{ } Since the atomistic results agree with LFM predictions that cover crack behavior also on macroscale, realization of the presented micro-mechanism is probable also in macroscopic specimens.  Figure 10 illustrates that the loading is sufficiently slow to achieve an agreement with the static LFM prediction at an applied stress level of 1 GPa without any damping in Newtonian equations of motion. In other words, dynamic effects in MD are negligible in this case and figure 10 confirms that the same boundary correction factors can be used for the atomistic 3D sample as predicted by LFM for macro-specimens. The other stress components have a similar character as in figure 3.
With increasing linear loading 1 PH, the crack was initiated in the middle of the crystal (layer 16) at the applied stress σ i =3.013 GPa. Corresponding stress intensity is   Figure 12(a) shows a plastic zone created in front of the crack in our observation plane (110), layer 16. Figure 12(b) is a detail in a perpendicular plane 110 (¯) and it clearly shows that the plastic zone is caused by twinning on the oblique planes (112). No significant cross slip of the emitted dislocations to {112} planes was monitored in MD at temperature of∼0 K and loading rate 1 PH. The oblique twins in front of the crack removes the LFM stress concentration at the crack, which hinders the further crack advance. The final character of fracture under slow loading 1 PH is more ductile, which illustrates the larger value of COD f ∼123 Å in comparison with the faster loading 5 PH, where COD f ∼40 Å.
To assess generation of the oblique twins in the oblique slip system 111 112 [¯]( ) by the crack itself we may use the LFM prediction (5) and consider a distance r=b=2.4825 Å and K I =K MD =0.871MPa m 1/2 corresponding to RK max mentioned above. We obtain larger slip stress τ b than the stress barrier τ twin for twin generation presented in [24] for the used potential: / since 3-layers twin is the minimum for its stability [26,27]. This short discussion shows that oblique twinning on planes {112} in front of the crack is possible also by continuum predictions, which implies its possibility also in experimental macro-specimens. Note that twins can arise also via dissociation of the core both in the case of the screw [28] and edge dislocations [29]-see e.g. Figures 9(a) and (c).
MD simulations at∼0 K also show that the threshold values for emission of the oblique dislocations á ñ 111 011 { }and oblique twins at the crack front are somewhat different: The threshold values can be changed by thermal activation in the system which is presented in the next section.

3.2.
MD results at the room temperature and their comparison with fracture experiments As mentioned in section 2, loading in this section starts after the previous heating at the time step 2000 h in all presented simulations at 295 K-see figure 13. In this section we will compare MD results for a longer crack a=64 d 110 under fast loading 5 PH and the short crack a=46 d 110 under slower loading 1 PH, which is relevant for fracture experiment in [3].  { } i.e. the same mechanism of crack growth was monitored both at 0 K and at 295 K. This is illustrated by figure 14 for time 13500 h. When the horizontal slip traces from screw dislocations á ñ 111 112 { }cross the crack plane (001) on inner layers (110) (e.g. the layers 6, 16, etc) they prepare a pre-processing zone for easy crack extension via separation of the planes (001) -see figure 14(b) where the position of the original crack tip point is denoted via the arrow. Cleavage crack extension in middle layer 16 at later time is illustrated in figure 15(a). However, generation of the oblique twinning on planes {112} was observed in front of the crack before the final fracture ( figure 15(b)), unlike corresponding simulations at∼0 K in section 3.1.1. When the twins reach the right free surface BL, new stress concentrators arise at free surface BL and new oblique dislocations á ñ 111 011 { }are emitted from the free surface. This increases crack opening and the final phase of fracture is realized via a small plastic necking of the atomistic sample. The final phase of fracture might be influenced by the finite sample dimensions. However, it did not occur at 0 K in figure 4 under the loading 5 PH and it may not occur even at 295 K, as mentioned below.
Oblique twins occurred in all runs 5 PH we tested at elevated temperature (see table 1), but randomly at different times, which affected the velocity and character of crack extension. The results are summarized in table 1. The largest portion of cleavage fracture was observed at the loading 1.59 GPa/10067 h and at 1.23 GPa/ 7788 h, which illustrates a 3D visualization of the fracture surfaces in figure 16. Table 1 is also used to assess thermal activation in the system. As to surface dislocation emission á ñ 111 011 ,    { }with the Burgers vector b/3 and thus, we may consider thermal activation of∼ k T 30 B for bulk crystal from [11]. Using K twin from relation (7), then  table 1). This caused that no fracture was monitored up to time step 50000 h. Note that a small change in loading rate to dP/dt=3 kN s −1 improved the situation and the overall shape of the atomistic sample at the final fracture was similar to the brittle fracture in figure 4 at 0 K. Nevertheless, the inspection of fracture surfaces showed that more cleavage facets occurred in MD simulations using the standard loading rate 2.92kN s −1 (5 PH), namely the ramp loading 1.23 GPa/7788 h-see figure 16. Table 1 shows also other differences between the simulations at 0 K and 295 K. As mentioned in section 3.1.1, the maximum of the resulting force RK in one half of the crystal physicaly represents a global resistence of the atomistic system against fracture and it occurs on the time scale in each case before the dislocation emission and crack initiation. Consequently, σ MD represents an internal critical stress available in MD at elevated temperature, K MD is global fracture toughness in MD at elevated temperature and G MD is the corresponding energy release rate in the heated atomistic sample. All the MD quantities are larger than σ i (or σ e ), K i and G i for crack initiation, derived from external loading, due to thermal activation. The average value of G MD at 295 K is 5.35 J m −2 and the difference -     { }was observed. Soon after crack initiation the oblique twins were generated-see figure 17. As a consequence, ductile fracture was monitored in MD, which illustrates figure 18 with the deformed fracture surface after the oblique twining. Such a sloping fracture surface is presented in figure 7 in [4], coming from SEM analysis of the fractured specimen No10 with the edge crack a/W∼0.3 loaded slowly in mode I.
Note that the oblique twinning and oblique dislocation emission á ñ 111 011 { }at the edge crack (001)[110] were already briefly reported in [32], where 3D simulation of the fracture experiment [4] on the shorter SEN specimens is presented. Oblique emission á ñ 111 011 { }was also reported in [23] under cyclic loading in mode I and mode II, where it resulted in crack deflections. No cross slip of the screw dislocations á ñ 111 011 { }and their contribution to crack extension is reported in [23,32]. Table 2 also includes comparison of the slow ramp loading 1 PH and the fast ramp loading 5 PH. The results presented in table 2 and the comparison of fracture surfaces in figure 18 for 1 PH and the fracture surfaces in figure 16 for 5 PH show that the slow loading 1 PH supports more oblique twinning, which lead to a ductile fracture. To verify that, we also performed MD simulations under linear loading in time, briefly presented in the next section.  [14]) from the loaded borders to the x-axis (where they are generated). In our case Δt x =1162 h, which is not negligible namely under faster loading 5 PH. The second disadvantage is that the high stress level before fracture makes the small 3D atomistic samples more ductile than expected in relation to experiments. Nevertheless, the initial crack behavior in the vicinity of the initiation can be described qualitatively well in relation to the experiments.
Comparison of the corrected MD results at temperature of 295 K with linear loading 1 PH and 5 PH is presented in table 3. Under linear loading the values COD f and W f /BW are larger than those in table 2 (ramp loading) because of the larger external forces acting during linear loading. This supports more oblique twining in the atomistic samples, which leads to more ductile behavior in comparison with ramp loading. Nevertheless, the values K i , COD f and W f /BW presented in table 3 show that the slow loading 1 PH leads to more ductile behavior in comparison with 5 PH, which is in a qualitative agreement with fracture experiments [3,4]. In this study it is illustrated in figures 19(a) and (b) for the linear loading 1 PH and 5 PH shortly after the crack initiation. Figure 20 shows the atomistic dependencies RK-t and RK-COD, in analogy with the experimental dependencies P-t and P-ΔL presented in [3,4] and in other experimental studies, reviewed e.g. in [2]. As already mentioned in Introduction, one may expect only a qualitative agreement between MD and the experimental results, not quantitative. And really, in agreement with the experimental results on bcc metals reviewed in [2], our MD curves RK-t and RK-COD with the increasing quasi-static loading rate are displaced to the left and towards the lower peak values, similar to experiments, i.e. with the increasing loading rate the fracture is faster and less ductile, or more brittle respectively. The same is valid for the ramp loading 1 PH and 5 PH that is more informative to answer our question on the contribution of the crack itself to ductile versus brittle behavior under different loading rates. The faster loading rate 5 PH leads to more brittle fracture, while the slower loading rate 1 PH to more ductile fracture, which was confirmed via inspection of the fracture surfaces in MD under the loading 1 PH (figures 18) and 5 PH ( figure 16). The presented MD results comply with our experiments [3,4] that motivated this study. The grape shape plastic zone in figure 4 from the experiment [3] and in figure 2 from [4] at the short crack a/W∼0.3 under the loading rate 1 mm min −1 is similar to MD figures 12(a) and 19(a) showing the oblique twinning in front the short crack a/W∼0.3 under the slow loading 1 PH. Under the fast loading 5 PH (for 5 mm min −1 ) such a plastic zone does arise neither in MD nor in the experiment. The average velocity of the crack growth under the ramp loading 5 PH is Vcr=417 m s −1 , as follows from table 2. If we suppose such a crack speed in the experiment [3] with the longer crack a/W=0.422 (where W∼10 mm and loading rate was 5 mm min −1 ), then the time needed for fracture will corresponds to about 1.4×10 -5 s, which is too short time for monitoring the crack extension via classical optical microscope and the process seems like a sudden fracture, as observed in [3].

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Funding-Acknowledgments
The work was supported by the Institute Research Program RVO: 61388998, further by the projects of the Technology Agency of the Czech Republic under grant TK01030108 and of the Czech Science Foundation with grant number GA19-03282S. The work of P. Hora was supported by the European Regional Development

Declaration of competing interest
We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.