Andrew F Calder
email: afcalder@liv.ac.uk

David John Bacon

Alexander V Barashev

Yuri N Osetsky

On the origin of large interstitial clusters in displacement cascades

Keywords: molecular dynamics, radiation damage user supplied): displacement cascades, self-interstitial atom clusters, shock waves radiation damage, displacement cascades, self-interstitial atom clusters, shock waves, molecular dynamics

Displacement cascades with wide ranges of primary knock-on atom (PKA) energy and mass in iron have been simulated by MD. New visualisation techniques are introduced to show how the shock-front dynamics and internal structure of a cascade develop over time. These reveal that the nature of the final damage is determined early on in the cascade process. We define a zone (termed 'spaghetti') in which atoms are moved to new lattice sites and show how it is created by a supersonic shock front expanding from the primary recoil event. A large cluster of self-interstitial atoms can form on the periphery of the spaghetti if a hypersonic recoil creates damage with a supersonic shock ahead of the main supersonic front.

When the two fronts meet, the main one injects atoms into the low-density core of the other: these become interstitial atoms during the rapid recovery of the surrounding crystal. The hypersonic recoil occurs in less than 0.1 ps after the primary recoil and the interstitial cluster is formed before the onset of the thermal spike phase of the cascade process. The corresponding number of vacancies is then formed in the spaghetti core as the crystal cools, i.e. at times one to two orders of magnitude longer. By using the spaghetti zone to define cascade volume, the energy density of a cascade is shown to be almost independent of the PKA mass. This throws into doubt the conventional energy-density interpretation of an increased defect yield with increasing PKA mass in ion irradiation.

Introduction

Exposure to fast neutron irradiation in the core of a nuclear reactor can change the dimensions and mechanical properties of materials. These effects are due to microstructural changes that evolve with time and arise from the production of radiation damage. The primary damage is produced in displacement cascades by the recoil of primary knock-on atoms (PKAs) and consists of point defects and their clusters. The average number of vacancies and selfinterstitial atoms (SIAs) per cascade and the fraction of them created in clusters vary with PKA energy and mass, temperature and material. In recent years, much research in the field of radiation damage has focused on SIA clusters, in recognition of the fact that their thermal stability, one-dimensional mobility and interaction with impurity atoms are key elements in explaining several long-standing problems, and may provide a strong driving force for the damage accumulation [1].

Although the presence of small clusters of SIAs in displacement cascades has long been inferred from experiments, e.g. [2,3], it has not been possible to determine the mechanisms involved by direct experimental observation because the relevant length and time scales are measured in nanometres and picoseconds. Our present understanding of the formation and properties of defect clusters in cascades has largely been provided by molecular dynamics (MD) computer simulation. MD simulation of model crystals of sufficient size to contain cascades became possible in the late 1980s [4][5][6][7][8] due to increasing computer power and development of many-body interatomic potentials [9,10], which provide realistic description of many important material properties yet offer computational efficiency. MD simulation has now become an invaluable source of information on cascades.

MD studies have shown that the number, N F , of vacancies or SIAs that survive intra-cascade recombination is fairly insensitive to metal and temperature. It increases with PKA recoil energy, E PKA , but, when expressed as a fraction of the NRT standard value [11], is approximately independent of E PKA for energy above about 5 keV [12,13]. It is also insensitive to the presence of solute atoms, such as substitutional Cu [14,15] and Cr [16] and interstitial C [17] in α-iron. This is because the number of defects produced in a cascade is governed by early stages of cascade development, when the atomic energies involved are much higher than the binding energy between atoms. The picture that has emerged regarding temperature and E PKA [13,[18][19][20]. Furthermore, despite the recognised importance of intracascade clustering for radiation damage effects in metals, understanding of the actual mechanism by which SIA clusters form and survive in a cascade has not been forthcoming. This is the issue addressed in this paper.

New insight into the cascade process has arisen as a by-product of another MD study of α-Fe [21,22], which in turn relates to experimental research undertaken on Fe and other metals some years ago, e.g. [23,24], and more recently [25,26]. In these experimental studies, selfion and heavy-ion irradiations of thin foils were used to investigate the effects of cascade parameters E PKA and PKA mass, m PKA , on damage. Transmission electron microscope (TEM) examination of the foils showed that an increase of the incident ion mass increased the number of visible clusters formed per ion and the mean loop size. (The clusters were vacancy dislocation loops: SIA defects were assumed to be lost to the surfaces.) Thus, although N F given by the NRT standard depends principally on E PKA and the displacement threshold energy of the metal, PKA mass clearly affects defect production and clustering in cascades.

In order to explain this, it was argued on the basis of binary-collision estimates obtained from [START_REF] Winterbon | Matematisk-Fysiske Meddelelser Udgivet af det[END_REF], that the cascade volume is decreased by an increase in m PKA , with the result that higher energy density and vacancy concentration are achieved for the same E PKA . This produces a higher probability of cascade collapse and vacancy loop formation.

The PKA mass effect in bulk Fe has been investigated recently by MD simulation of cascades created by recoil of atoms with mass varying from 12 amu (carbon) to 209 amu (bismuth) and energy in the range 5-20 keV [21]. It was found that large vacancy and SIA loops are produced with increasing frequency as PKA mass increases. These studies have now been extended to other conditions, and the structure and statistics of the cascade damage are the subject of another paper [22]. However, by using new visualisation techniques to probe the stages at which damage is formed in these high m PKA cascades, it became possible to reveal the mechanism by which large clusters of SIAs actually form early on in the cascade process. Furthermore, a meaningful definition of 'cascade volume' has emerged. The insight gained into the cascade process is another topic of this paper. The conclusions are drawn in section 7.

MD simulations

The MD model of Fe was the same as that described in [21], but with a wider range of cascade parameters. In summary, it contained between 0.25 and 2 million atoms, depending on E PKA , which ranged from 5 to 30 keV. Atomic interactions were computed using the embedded-atom-method empirical potential from [28], which was developed with specific emphasis on correct description of the relative stability of SIAs of different configurations and the vacancy migration energy. The short-range (≤ 0.1 nm) pair part of the potential was replaced by the Universal screened-Coulomb potential [29] in order to better describe the high-energy atomic collisions that occur in cascades. The two potentials were joined smoothly by an exponential spline over the range 0.100-0.205 nm. Electron-atom interactions were not included.

The crystal was equilibrated at 600 K prior to initiation of a cascade and maintained at constant volume with the lattice parameter, a 0 , corresponding to the zero-pressure condition.

Integration of the equations of motion of atoms was performed using a leapfrog algorithm with variable time-step, determined by limiting the displacement of the fastest atom to 0.003a 0 in one time-step. No temperature control was applied and the temperature increase due to introduction of the PKA was less than 70 K. Fe, Bi and 'molecular' Bi 2 and Bi 3 . For C and Bi, one Fe atom was substituted. For the 'molecular' PKAs, two (or three) second nearest-neighbour atoms were substituted and used as PKAs, each with the same direction and half (or one third) of E PKA . At least 30 cascades were simulated for each combination of PKA energy and mass. The interatomic potentials for the C-Fe and Bi-Fe interactions were the same as that for Fe-Fe except at short range, where the Universal screened-Coulomb potential for the required atomic pair was used for the pair part, as described in the preceding paragraph. We note, however, that the PKA-Fe shortrange potential has little effect on cascade damage compared with PKA mass [21].

Cascades in bulk

Spaghetti analysis and visualisation

Unambiguous identification and cataloguing of point defects and their clusters is challenging when a large number of cascades has to be treated. Effort has tended to concentrate on two methods to identify damage. One is based on analysis of Wigner-Seitz (W-S) cells and is used to identify vacant sites (≡ cells with no atoms) and SIAs (≡ cells with more than one atom). This method is adequate for determining the number, N F , of surviving vacancies (equal to the number of surviving SIAs), but is not appropriate for identification of structures with correlated atomic displacements typical of stacking faults and interstitial dislocation loops with perfect Burgers vector [30]. The second method is displaced-atom (D-A) analysis, based on the position of atoms with respect to lattice sites. For each lattice site the closest distance to an atom is evaluated and if larger than a threshold value (0.3a 0 here) the site is treated as vacant. Similarly, for each atom the closest distance to a lattice site is determined and, if larger than the threshold value, the atom is treated as 'displaced'. Evaluation of the number of point defects is not straightforward in this approach, either in the ballistic phase of a cascade, when many atoms undergo temporary elastic displacements, or at the end, when a dumbbell SIA is flagged as a vacant site and two displaced atoms, and a crowdion is revealed as a line of displaced atoms interspersed with one fewer vacant sites. Both methods were applied in the present work.

Figure 1 shows the variation of the number of displaced atoms calculated by the W-S cell and D-A methods versus time for a 20 keV Bi 3 cascade, for which N F = 77. The number given by the D-A analysis is an order of magnitude larger at all times, for the reasons given above. The To complement these two treatments of the vacancy and interstitial components of the primary cascade damage, we have also analysed the final coordinates of atoms in order to identify those that occupy lattice sites other than their original ones. To visualise the movement of atoms in this category, it has been found convenient to construct a line between the original and final sites of each atom. In a cascade involving permanent displacement of hundreds or thousands of atoms, nearly all site-to-site lines form a chain of connectors between sites in a sequence of atomic replacements, and there may be many such chains in the cascade zone. It should be noted that a chain is not the trajectory of a single atom. It is simply the path of many atomic replacements leading from a vacancy to a final SIA site. For want of a better word, we call such a set of chains 'spaghetti'. Some possibilities are illustrated schematically as follows.

Figure 2(a) shows a single chain starting at vacant site 'A' and ending at atom 'E', which is part of the single dumbbell SIA consisting of atoms 'E' and 'G' and vacant lattice site 'F'.

The atom originally at site 'A' was displaced to the site 'B', the atom originally at 'B' was displaced in turn to site 'C', and so on. fairly typical of those that arise from PKAs of this type and energy, but the number of defects was the smallest among the 30 events simulated, namely N F = 4 compared with a mean of 37.8 and a maximum of 82. The damage consists of four single vacancies, one single SIA and a 3-SIA cluster. Thus, despite the apparent complexity of the spaghetti -it has 8856 site-tosite links in the chains -the number of chains involved is only eight, as drawn schematically in Figure 3(c), where the interstitial cluster is shown in the shaded region in the centre of the Figure . Visualisations of all eight chains are presented in Figure 4(a-h).

[Insert Figure 3 The principal features that stand out in these and similar Figures and their analysis are as follows.

[Insert Figures 5, 6 and 7 about here] (i) The spaghetti zones in the light-mass C-PKA cascades (Figure 5) are diffuse and usually dispersed into sub-zones. These events with small m PKA produce a significant proportion of their interstitials as single SIAs outside the spaghetti regions. An atom in one of these defects is joined to its parent spaghetti by short links associated with migration of the SIA, but the continuation of the chain to a vacant site within the spaghetti may be complex.

(ii) The spaghetti produced by the high-mass Bi 2 PKAs (Figure 7) is denser and rarely exhibits sub-zones. The number of single SIAs linked by collision and migration chains to the main zone is relatively small for these cascades. The PKA event is always at the edge of the zone. Vacancies are created within the spaghetti and SIAs are formed at its periphery.

(iii) The self-ion cascades (Figure 6) tend to fall between the extremes of (i) and (ii). Fewer spaghetti sub-zones are formed than in the C-PKA cascades, but the probability of a single-zone cascade producing a sizeable SIA cluster is less than in the high-m PKA damage. (vi) The probability of formation of clusters of interstitial atoms increases with increasing PKA mass. (The nature of these clusters and the vacancy arrangements within the spaghetti zone core is the topic of another paper [22].) However, the creation of a dense spaghetti does not guarantee formation of a cluster. Thus, the formation of defects requires conditions other than simple deposition of energy, which is revealed in the huge mixing shown by the spaghetti. Clustering of SIAs increases defect survival and the conditions leading to clustering are discussed in the next section.

Cascade dynamics

With reference to the plot in Figure 1, the rapidly expanding zone of displaced atoms exhibits two distinct stages before the number of displaced atoms reaches the peak at time t ~ 1 ps. Up to t ~ 0.1 ps, the number of displaced atoms in the D-A treatment is approximately proportional to t 3 (dotted line), implying a zone front moving with approximately constant velocity. Estimates obtained from either visualisations of displaced atoms or the change in volume of the zone with time show that the front moves outwards with a speed of approximately 30 km s -1 , i.e. about five times the longitudinal velocity of sound, v s , in iron.

With the elastic constants calculated for the potential used here, v s varies 5.6 to 6.5 km s -1 at 0K for different low index directions. The variation is 5.1 to 6.1 km s -1 with the experimental elastic constants measured at 595K [31]. Thus, the front is supersonic in character, i.e. speed in the range v s to 5v s , up to t ~ 0.1 ps, but becomes sonic as t increases further to 1 ps. More information on these supersonic and sonic shock fronts follows later in this sub-section.

In the context of the present paper, W-S visualisation of interstitial atoms and vacancies of some cascades reveals features that are characteristic of interstitial cluster formation. Figure 8 contains images generated by the W-S visualisation method of a 30 keV Bi 2 cascade at times from 0.04 to 32.6 ps, as indicated. The image in Figure 8(a) for t = 0.0402 ps shows that two hypersonic recoil (HR) events labelled HR1 and HR2 emerge ahead of the main zone of damage. Their speed is near to 10v s , i.e. they are hypersonic in character (speed > 5v s ). These hypersonic events start to create their own damage zones, which develop further and generate supersonic shock waves ahead of the primary shock front by t = 0.1 ps (Figures 8(b) and(c)).

The shock fronts meet and the damage zones impinge on each other by about 0.2 ps (Figures 8(d) and (e)). The spaghetti zone and the W-S images are superimposed in Figure 8(e) (t = 0.2558 ps) Figure 8(f) shows that by the end of the cascade process, a large interstitial cluster is formed in the region where the hypersonic recoil damage was created. Analysis of all the cascades simulated in the present work shows that subsidiary damage due to hypersonic recoils produced early on (t < 0.1 ps) in the ballistic phase is associated with the existence of large clusters of SIAs in the final damage state. This statement is examined further in section 5, but before that we consider the nature of the shock processes in more detail.

[Insert Figure 8 about here] The terminology shock front implies an abrupt change in pressure, density and temperature.

We have examined the latter two parameters in the following way. Figure 9 The sequence of plots in Figure 9 show that, for t up to about 0.1 ps, T has an almost discontinuous change at the supersonic front between the cascade core and the surrounding crystal. Behind the front, crystalline order is destroyed and atoms have temperature in excess of 2000 K, which is just above the melting point of 1750±25 K for the interatomic potential used [32]. This hot core zone also has regions near the centre with no atoms. The shock front becomes less abrupt and its velocity slows by t = 0.2 ps (Figure 9(c)), beyond which there is noticeable diffusion of thermal energy into the surrounding crystal. The 'molten' zone reaches its maximum size by 0.3 to 0.4 ps, which is much shorter than the time of the peak in the number of displaced atoms in Figure 1, but it only starts to shrink after 3 to 4 ps. The region of spaghetti defined in section 3 corresponds to this zone at its maximum size. This description of a highly disordered, hot cascade core is consistent with the cross-sections showing atomic positions produced from early computer simulations of displacement cascades, e.g. [4][5][6].

Another way to visualise the shock-wave nature of the disturbance created by a cascade is to consider the atomic density. This is achieved in Figure 10 by plots showing the atomic neighbour density in the same cross-section and at the same times as those in Figure 9. The neighbour density of each atoms is defined to be the number of atoms in a sphere centred on the atom in question. A sphere radius of 1.405a 0 was found to be suitable because it excludes the twelve third-nearest neighbour atoms at a distance of 1.414a 0 in a prefect crystal, but is sufficiently close to it to enable increases in density associated with shock effects to be detected. The neighbour density is seen to be high (exceeding 26 at some sites) at the shock front, which is compensated by very low values (between 0 and 8) in the central core of the cascade. The front defined by density is more obvious than that in Figure 9 for time greater than about 0.2 ps, when it has sonic speed and the deformation associated with it is purely elastic.

[Insert Figure 10 about here] It was noted above that W-S identification of displaced atoms reveals two hypersonic collision events created ahead of the supersonic front in this cascade. They are apparent as two 'ears' at about one o'clock in Figures 9(a) and (b), and appear to be subsumed within the zone behind the main front in Figure 9(c). However, Figure 10 shows that they actually retain a separate identity and create their own shock fronts before being reached by the main front

On the origin of large SIA clusters

As explained in section 1, the computer simulations in [21] showed that the probability of producing large vacancy and SIA clusters in a displacement cascade increases with increasing PKA mass and energy. This study has now been extended and results on the structure and statistics of the clusters are the subject of another paper [22]. The relevance of the present work is that by using the visualisation techniques described to examine cascade development and the final damage, it has been possible to observe a correlation between the production of large clusters of SIAs and a process initiated at a very early stage of cascade development.

Correspondingly, the creation of large SIA clusters by events early on in the cascade process leads to the formation of large vacancy clusters by the end of the thermal spike phase, i.e. after many picoseconds.

We have observed that the formation of an SIA cluster that becomes distinguishable as such by the end of the collision phase, i.e. at times of a few picoseconds, proceeds by creation of a hypersonic energetic recoil atom that travels beyond the main shock front and creates its own shock wave and associated damage near the edge of what later becomes the spaghetti zone. In other words, the site for the cluster is created in less than 0.1 ps from the time of the initial primary recoil. (The time when the hypersonic recoil stops is indicated in Figure 1). The process is summarised schematically in Figures 11(a)-(d). Most of the supersonic shock front does not encounter damage formed by a hypersonic event and the atoms displaced by the front outside the spaghetti zone simply return to their lattice sites by elastic relaxation after the time of the peak in Figure 1. However, when an energetic impact event produces one or more hypersonic recoils (two are shown in Figure 11(a)) beyond the main front (see Figure 11(b)), the main and secondary supersonic shock waves with high atomic density at their fronts and low density disorder behind expand and meet (Figure 11(c)). It appears from the evidence of the visualisation techniques that the main front then injects atoms into the core of the secondary region(s) (Figure 11(d)). The number of atoms there exceeds the number of lattice sites and the excess atoms become trapped by the rapid recovery of the surrounding crystal (Figure 11(e)) and form an interstitial cluster. (Note that it seems unlikely that interstitial atoms form in the compression region where the two fronts meet, because the Burgers vectors of the interstitial dislocation loops that eventually form are found to be tangential to the spaghetti zone [22], as depicted in Figure 11(d).) The interstitial atoms formed in this way are matched by vacancies formed during solidification of the hot disordered core of the main zone, which is a relatively slow process. It was inferred from atomic configuration crosssections produced in early simulations of cascades in fcc metals that free volume accumulates on the liquid side of the solidification interface, with the result that a vacancy cluster forms in the centre of the cascade zone, e.g. [4][5][6].

[Insert Figure 11 about here] It should be noted that the formation of dislocation loops seen here is not the same as the prismatic loop punching mechanism postulated by Diaz de la Rubia and Guinan [33], which would produce loops with Burgers vector orientated radially from the cascade core. Nor does it appear to be related to the mechanism of shear dislocation loop creation behind planar shock fronts produced by impact. Germann et al. [34] and Kadau et al. [35] have used MD to simulate the effects of a planar shock front in models of fcc and bcc crystals, and interpreted the dislocation structure observed in terms of earlier models of shock effects [36][37][START_REF] Armstrong | ONR Workshop on Dynamic Deformation, Fracture, and Transient Combustion[END_REF]. The dislocations were created in the strained crystalline material behind the front, which is not the process found in the present work. In the simulations of iron crystals, Kadau et al. [35] found that solid-solid, martensitic-like transformations were induced by shock. Transformations were not present in the cascade damage regions here, presumably because the melt-like zone has a long enough life and is sufficiently small for it to solidify coherently with the surrounding bcc crystal. Thus, we have established a connection between formation of large SIA clusters and the existence of a high-energy, hypersonic recoil events in the primary stage of cascades. We first became aware of the connection through visualising effects in cascades due to PKAs with high mass, and it then became apparent in all the cascade simulations that resulted in large clusters. It seems essential that the hypersonic recoil atom stops beyond, but near, the edge of the spaghetti region. If it stops well beyond the supersonic front, it results in sub-cascade that develops independently of the main cascade. If it were within the vacancy-rich cascade core, where the energy deposited is much higher and the lifetime of such an excitation is long, the damage would be effectively annealed.

On energy density arguments in cascade collapse

It was noted in section 1 that the effects of parameters E PKA and m PKA on cascade damage on Fe and other metals, e.g. [23][24][25][26], have been investigated experimentally by ion irradiation of thin foils and subsequent examination of the visible defects, i.e. dislocation loops, in a TEM. This has shown that increasing ion mass increases the number of visible clusters per ion and the mean cluster size. Thus, although N F given by the NRT standard depends principally on E PKA , PKA mass clearly affects defect production and clustering in cascades. In order to explain this, it is generally argued on the basis of binary-collision estimates that the cascade volume is decreased by an increase in m PKA , with the result that higher energy density and vacancy concentration are achieved for the same E PKA . This produces a higher probability of cascade collapse. Note that collapse creates vacancy loops late on in the cascade process, not interstitial loops. SIA defects are assumed to be lost to the nearby surfaces in these thin-foil experiments.

The present work offers a different explanation for the increasing probability of formation of visible vacancy loops with increasing PKA mass. The supersaturation of vacancies in the cascade core, which is necessary for cascade collapse to a vacancy loop, has been shown to arise from the prior clustering of interstitial atoms, which is initiated by a hypersonic recoil process in less than 0.1 ps from the time of the primary recoil event. Furthermore, the techniques of visualisation presented here allow the concept of 'cascade volume' and 'energy density' to be considered from a more realistic standpoint than that adopted hitherto. The spaghetti analysis permits us to define an effective region within which permanent damage occurs, and then calculate its volume and, hence, energy density. Calculation of the volume was performed using the MSMS code of Sanner, Spehner and Olson [START_REF] Sanner | [END_REF], which is normally used to calculate molecular surfaces and the solvent-accessible surface (SAS) and solvent-excluded surface (SES) of large protein and nucleic acid molecules. For the cascade calculations it was assumed that the spaghetti region consists of a set of spheres representing iron atoms centred on the 'joints' in the spaghetti chains with radius r s and a MSMS probe sphere radius of radius r p , as shown schematically in Figure 12(a). The SES topological boundary was chosen as the appropriate surface of the spaghetti region surface area and associated enclosed volume. Values r s = 0.55a 0 , i.e. slightly larger than the hard-sphere atomic radius of 0.42a 0 , and r p = 0.6a 0 were found to be suitable, although some surfaces required a slight adjustment to r p to allow the SES to be computed successfully. If spatially-separated SES volumes were identified for a single cascade event, the sum of all the identified SES areas and volumes was taken. Visualisations of the MSMS-computed SES were done using PyMOL [START_REF] De Lano | The PyMOL Molecular Graphics System[END_REF] and a Python script msms_cgo.py from [START_REF] Campbell | Queen's University[END_REF]. (PyMOL was also used for the D-A, W-S and spaghetti images shown in this paper.) An example for the cascade of Figures 8-10 is presented in Figure 12(b), which is a superposition of the interstitial atoms and vacancies using W-S visualization, the spaghetti and the MSMS probe spheres. It can be seen that the volume and surface area include the main body of spaghetti as well as chains that lead to single SIAs (e.g. top right of figure) and simple replacement sequences that emerge from, and return to, the main body. For cascades with high m PKA , the contribution of such chains is small.

[Insert Figure 12 about here] Figure 13 contains mean values (with bars representing the standard deviation) for the area, A s , and volume, V s , of the spaghetti zone of all the cascades for each of the combinations of E PKA and m PKA considered here (619 cascades in total). Several features stand out. First, V s depends approximately linearly on the PKA energy. Second, the surface area of the region of permanent damage is seen to be sensitive to the PKA mass and decreases with increasing mass. This reflects the more irregular shape when m PKA is reduced. Third, and most important, the cascade volume, and hence energy density, is almost independent of PKA mass. This is in stark contrast to the conventional understanding derived from binary collision models, and leads to the conclusion that interpretation of trends from experimental measurements of vacancy defect yield based on energy density arguments is not applicable for the analysis of cascade collapse.

[Insert Figure 13 about here] 3. In some cascades, a large cluster of self-interstitial atoms is formed on the periphery of the spaghetti. This is the result of interaction between two kinetic elements, namely a hypersonic recoil that occurs in less than 0.1 ps after the primary recoil and creates damage with its attendant supersonic shock wave ahead of the main supersonic front.

When the two shock waves meet, the main one injects atoms into the low-density core of the other. The number of atoms there then exceeds the number of lattice sites available during the rapid cooling and elastic recovery of the surrounding crystal.

4. Large interstitial clusters are therefore formed before the onset of the thermal spike phase, i.e. in less than a few tenths of a picosecond after the primary recoil. A correspondingly large number of vacancies forms later in the core of the spaghetti as the crystal lattice is restored during cooling, i.e. at times one to two orders of magnitude longer.

5. By defining cascade volume in terms of the permanently damaged spaghetti region, the energy density of a cascade is shown to be independent of the PKA mass. This result rules out the conventional energy-density interpretation of the increased probability of defect cluster formation with increasing PKA mass observed after ion irradiation.

3 SIA

 3 clustering is less clear. The clustered fraction of SIAs and the structure of the clusters formed by the end of the cooling down phase of the cascade region have been shown by MD to depend on many factors, including crystal lattice structure, stacking fault energy,

  MD simulations are first summarised in section 2 and then the results of a novel technique developed to visualise the displacement of atoms in a cascade are described with examples in section 3. As explained in section 4, the analysis assists understanding of the link between the dynamic processes associated with shock fronts created in the collision phase of cascades and the hot zone that characterises the thermal spike phase. This leads in section 5 to appreciation of the origin of large SIA clusters formed in cascades created by PKAs with high energy and mass. It is demonstrated in section 6 that the conventional energy-density arguments used to explain the formation of vacancy loops by cascade collapse are invalid.

  peaks due to shock waves re-entering the MD cell through the periodic boundaries, although the W-S plot shows that these have no effect on the final damage. The regions (with labels) identified in Figure1are discussed in section 4.[Insert Figure1about here]

  Figure 2(b) shows two such chains, each of which gives rise to a Frenkel pair. Figure2(c) presents a possibility in which three chains produce only two Frenkel pairs. Chains 2 and 3 both start with a vacancy and end at the same SIA, and chain 3 ends at the second SIA. Closed atomic replacement sequences, in which every atom replaces another and no defects are created, also occur.[Insert Figure2about here] The actual spaghetti zone and final damage arrangement produced by a 20 keV Bi 2 cascade are presented in Figures3(a) and (b), respectively. (The vacant sites and displaced atoms in these Figures are visualised by the D-A treatment. Most of the defects are obscured by the spaghetti in (a).) The site and direction of the PKA are indicated by the head and shaft of the arrow shown at the periphery of the zone at about the 10 o'clock position. This spaghetti is

  about here] [Insert Figure 4 about here] Samples of spaghetti and interstitial and vacancy damage (D-A visualisation) for 20 keV cascades produced by C, Fe and Bi 2 PKAs are presented in Figures 5, 6 and 7, respectively.

  The pure atomic mixing chains, i.e. those with no SIAs and vacancies, constitute only a small fraction of all displaced atoms. Estimates from several cascades suggest less than 10%.(v) It is clear from the N F data included in Figures 4 to 7 that a large amount of spaghetti is not necessarily associated with a large number of Frenkel pairs.

  contains snapshot cross-sections through the central region of the cascade in Figure 8 at various times from 0.04 to 3.50 ps. The cross-section plane is (01 1) and the direction of the abscissa and ordinate is [100] and [011], respectively. The scale unit of length is a 0 . The colour (grey shading in the printed version of the paper) indicates the temperature associated with each atom lying within 0.7a 0 of the plane of the Figure, as shown in the scale in Figure 9(a). Foran atom with kinetic energy E k , the temperature, T, is given by T = 2E k /3k B , where k B is the Boltzmann constant. To smooth extreme variations in the temperature from atom to atom, the temperature in the Figures is the average over all atoms in a sphere centred on the atom in question. A sphere radius of 1.1a 0 was found to provide reasonable sensitivity in the temperature plots. In a perfect crystal, a sphere of this radius contains the central atom and its 14 first-and second-nearest neighbours.[Insert Figure9about here]

11 (

 11 parts (c) and (d)). This is important for the formation of large clusters of SIAs, as shown in the next section.

1 .

 1 The new visualisation techniques introduced here show how the shock-front dynamics and internal structure of a displacement develop over time and that the nature of the final damage is determined early on in the ballistic phase.

2 .

 2 The spaghetti zone in which atoms are moved to new lattice sites by the cascade process is identified with a hot, disordered zone created as a supersonic shock front. The boundary of the spaghetti coincides with the transition region where the speed of the front slows to the sonic regime.

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2. Schematic representation of different types of chains of the spaghetti.

Figure 3 .

 3 Figure 3. (a) The actual spaghetti zone, (b) final damage arrangement produced by a 20 keV Bi 2 cascade and (c) schematic representation of chains. The vacant sites and displaced atoms are visualised by the D-A treatment.

Figure 4 .

 4 Figure 4. Visualisations of all eight chains of the spaghetti presented in Figure 3(c).

Figure 5 .

 5 Figure 5. Sample of spaghetti and interstitial and vacancy damage (D-A visualisation) for 20 keV cascades produced by C PKA.

Figure 6 .

 6 Figure 6. Same as in Figure 5 but for Fe PKA.

Figure 7 .

 7 Figure 7. Same as in Figure 5 but for Bi 2 PKA.

Figure 8 .

 8 Figure 8. W-S visualisation snapshots at various times of a cascade created by a 30 keV Bi 2 PKA. HR1 and HR2 in (a) are hypersonic recoils that emerge ahead of the primary shock front and result in the formation of the large SIA cluster indicated in (f). The spaghetti zone and the W-S image are superimposed in (e).

Figure 9 .

 9 Figure 9. Snapshot cross-sections through the cascades in Figure 8 at various times. The cross-section plane is (01 1) and the direction of the abscissa and ordinate is [100] and [011], respectively. The scale unit of length is a 0 . The colour (grey shading in the printed version of the paper) indicates the temperature associated with each atom lying within 0.7a 0 of the plane of the paper, as shown in the scale in the part (a). The temperature is the average over atoms in a sphere of radius r c = 1.1a 0 centred on the atom in question.

Figure 10 .

 10 Figure 10. The atomic neighbour density (the number of atoms in a sphere of radius r c = 1.405a 0 ) centred on the atom in question in the same cross-section and at the same times as those in Figure9.

Figure 11 .

 11 Figure 11. Schematic representation of cascade development over time and interstitial cluster formation.

Figure 12 .

 12 Figure 12. (a) Schematic illustration of the MSMS construction used to calculate the surface area and volume associated with a spaghetti chain. (b) Superposition of the interstitial atoms and vacancies using W-S visualization, the spaghetti and the MSMS probe spheres for the cascade shown in Figures 8-10.

19 Figure 13 .

 1913 Figure 13. Mean values with bars representing the standard error for the area, A s , and volume, V s , of the spaghetti zone of all the cascades for each of the combinations of E PKA and m PKA considered here.

  Fe were simulated by giving one or more atoms velocity (2E PKA /m PKA ) ½ in either <123> or <135> crystallographic directions. Five different PKAs were used, namely C,

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	5

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgments

The research was supported by grant GR/S81162/01 from the UK Engineering and Physical Sciences Research Council, and partly by the Division of Materials Sciences and Engineering and the Office of Fusion Energy Sciences, U.S. Department of Energy, under contract DE-AC05-00OR22725 with UT-Battelle, LLC.


Ghiath ; Monnet
email: ghiath.monnet@edf.fr

Yuri Osetsky

David John

Edf

David Bacon

Mesoscale thermodynamic analysis of atomic-scale dislocation-obstacle interactions simulated by molecular dynamics

Keywords: atomistic simulation, dislocation dynamics, dislocation interactions, dislocation theory, irradiation effects, thermal activation of deformation atomistic simulation, dislocation dynamics, dislocation interactions

Mesoscale thermodynamic analysis of atomic-scale dislocation-obstacle interactions simulated by molecular dynamics

Introduction

In a recent work [1], data on the mechanics of dislocation-defect interaction at 0 K obtained by molecular statics (MS) simulations in [2] were analysed in a mesoscopic framework. It was shown that the change in potential energy as a dislocation encounters and finally cuts through a periodic row of voids can always be expressed as the sum of three different components: (i) the energy of elastic deformation, (ii) the energy associated with the elongation of the dislocation as it bows out between obstacles and (iii) the interaction energy between the dislocation and an individual void. At finite temperature, the response of the system under mechanical load is no longer deterministic, for the local features of the interaction between a dislocation and obstacle introduce local energy fluctuations, even in isolated systems [START_REF] Reif | Fundamentals of Statistical and Thermal Physics[END_REF][START_REF] Kocks | [END_REF]. These fluctuations assist dislocations to unpin, thereby inducing an irreversible evolution of the whole system. This explains the general decrease of the obstacle strength with increasing temperature.

Although a large number of studies by molecular dynamics (MD) simulations of dislocation-obstacle effects have been reported in the literature, see for example [START_REF] Bacon | Dislocations in Solids[END_REF][START_REF] Marian | [END_REF]7], little has been done to analyse the nature of the localised thermal activation process. For example, recent investigations [8,9], although applied to the Peierls barrier rather than the localized obstacle case, provided two different treatments. Thus, analysis of the motion of a screw dislocation in iron in [8] based on a time-independent rate equation showed that the stress dependence on temperature can be correlated with the time-independent integration of the jump probability. This method provides a self-consistent definition of the critical stress that accounts for its dependency on strain rate and temperature. The properties deduced from this integration provided an explanation of the shift between the critical stress measured in MD simulation and experiment. In the MD simulation of the motion of the Lomer dislocation in aluminium in [9], another approximation is used to deduce the critical stress. The incompatibility of these two approaches is the subject of a forthcoming paper.

In a completely different approach, Jarzynsky [10,11] and Crooks [12] proposed a connection between equilibrium statistical mechanics and systems driven far from equilibrium within conditions of microscopic reversibility. The established theory, however, cannot be usefully applied in the case of thermal barriers.

It should also be noted that interpretation of the critical stress reached in MD simulations is not straightforward. For example, the connection with theory via the line tension approximation is not successful [13] because of the absence of an appropriate of strength defined by a critical cusp angle on the dislocation line, see for example [START_REF] Friedel | Dislocations[END_REF], cannot be used even if the local dislocation character is taken into account (see below) [13].

In this paper, we propose a new approach to analyse dislocation-obstacle interactions at finite temperature. We illustrate it for the case of edge dislocation-void interaction, as in [13]. The paper is organized as follows. First, we describe briefly in section 2 the features of the atomic-scale simulations used here. Then the analysis reported in [1] is applied and extended in section 3 in order to interpret the results of MD simulations. In this section we provide a new definition of the absolute strength of finite-size obstacles that is appropriate for mesoscopic treatment. The MD results of the dislocation-void interaction are presented and the effect of the temperature is discussed in section 4. The thermodynamic treatment of the interaction process is presented in section 5.

MS and MD methods for simulating dislocation-obstacle interactions

The source of the mesoscale analysis introduced in this paper is data from atom-scale modelling of a gliding edge dislocation in α-iron that interacts with an obstacle in the form of a spherical void centred on the glide plane. The edge dislocation with the Burgers vector b = ½[111] parallel to the x axis is initially aligned along the ] 2 11 [ direction of the y axis and glides on the ) 0 1 1 ( plane. We have used the two simulation techniques MS and MD with the interatomic potential developed by Ackland et al. [START_REF] Ackland | [END_REF], for which the empirical parameters were fitted for BCC iron. The atomic model contained 877,590 lattice sites and its x, y and z dimensions were approximately 25 nm (~100b), 21 nm and 20 nm, respectively. The void consisted of 59 vacancies and had diameter, D, equal to 1 nm. Details of techniques to create the dislocation, visualize its core atoms, induce its motion and calculate stress and strain are presented in [2]. Here we just mention a few details.

In the MS simulations, the dislocation moved due to resolved shear strain, γ, applied in increments ∆γ. The potential energy of the model was minimized by atomic relaxation after each increment, and the potential energy and the applied shear stress, τ, corresponding to the strain were then calculated for that set of atomic coordinates. The atomic positions obtained in this way represent the equilibrium configuration of the dislocation interacting with the void at each value of τ and γ, and are therefore the atomic-level equivalent to the dislocation line simulations were performed with ∆γ equal to either 10 -4 or 5x10 -5 .

Before proceeding, we recall the energy contributions analysed in the original paper [1]. In MS simulations, the incremental applied work δW app and the change in potential energy dU pot are coupled through the equation δW app = ∆U pot + δQ, where δQ is the heat extracted from the system in order to keep it at temperature T = 0 K. It accounts for all irreversible processes in the simulation [12], namely work against the lattice friction W pl and the annihilation of a dislocation segment inside the void. The interaction energy (E int ) can be deduced from the relation

E int = U pot -E curv -E el
, where E curv is the energy necessary to bend the dislocation before it unpins from the void and E el is the energy of the elastic deformation of the model due to the applied shear stress.

MD simulations were carried out for the temperature range 1 to 600 K at values of the applied strain rate, & γ , equal to 1x10 6 s -1 . Additional simulations for 600 K were also carried out for & γ = 5x10 6 s -1 . At these strain rates the steady state dislocation velocity is 2.0 and 10.0 ms -1 for the model size used. The total simulation time varied from 1.6 to 16 ns. During this time, the temperature of the model, calculated from the average kinetic energy of all the atoms, rose by less than 1 K. This increase consists of two contributions. One is due to friction in dislocation motion and depends on the dislocation velocity. The dislocation accelerates quickly when attracted into the void and when leaving it at the maximum stress.

The other contribution arises from the decrease in potential energy when the dislocation enters the void, losing part of its length, and leaves it, releasing the accumulated elastic energy. The latter is the biggest contribution in the cases considered but occurs after release of the dislocation and does not affect the dislocation-void interaction analysis. The local temperature increase due to potential energy change is an unavoidable yet realistic phenomenon. This thermal energy diffuses over the surrounding atoms but does not homogenize completely within the simulation time. There was practically no change in the temperature near the crystallite boundaries in these simulations and it was unnecessary to apply a thermostat.

For statistical treatment of the results, some simulation conditions were repeated up to five times but with different initial atomic velocity distributions. The principal outputs of the MD simulation for a given temperature are the dependence of applied stress, potential energy and kinetic energy on applied strain and/or time. Some results of an MS simulation (T = 0 K) are seen in Figure 1, where τ versus γ is plotted in (a) and the potential and other energy changes are plotted against γ in (b). The shear strain γ r denotes the reference state in which the dislocation is straight and penetrates the void through its centre. It corresponds to τ equal to the friction stress, τ f , for the edge dislocation, which is close to 23 MPa.

Fig. 1 : about here At this stress level, the potential energy is not minimum, for although the curvature energy is zero, the elastic energy is Vτ f , 2 /2µ, where V is the volume of the simulation box and µ the shear modulus. At γ r , E int , is equal to -7.5 eV, compared to the value -17.4 eV found in the data for a 2 nm void [1]. The former value is approximately half of the latter, demonstrating that the interaction energy, deduced from the analysis in [1], is a characteristic of the dislocation-obstacle interaction and independent of the simulations conditions. This interaction energy is basically needed to form the two steps on the void surface created when the dislocation enters and shears the void. However, the difference between -7.5 and half of -17.4 is due to the fact that the energy of the surface step on the 1 nm void is not simply half that of the step on the 2 nm void.

Critical strength of obstacles of finite size

It is clear that during pinning of the dislocation by the void, mechanical equilibrium prevails around the obstacle. Unpinning occurs when the obstacle experiences an effective force due to the dislocation curvature that is larger than it can sustain. The conventional picture of this equilibrium is that of point obstacles [START_REF] Friedel | Dislocations[END_REF][START_REF] Hirth | Theory of Dislocations[END_REF], depicted in Figure 2a. The line tension (LT), Γ, on each side of the obstacle generates a resultant force F = 2Γ cos(ϕ/2), where cusp angle ϕ reflects the obstacle strength. With an increasing line curvature, ϕ decreases and F increases. When the force exceeds the resistance of the obstacle, F max , the dislocation unpins from the obstacle.

Fig. 2 : about here This scheme has been used in theories of alloy strengthening, for a review see [START_REF] Nembach | Particle Strengthening of Metals and Alloys[END_REF]. It was also used for the treatment of the Orowan mechanism [START_REF] Bacon | [END_REF] by employing a line tension that accounts properly for the dipole interaction between dislocation arms around the particle.

However, the line tension approximation suffers from three principal limitations. (i) It is not easy to extend it to obstacles of finite size. The critical stress it gives with the actual ϕ for a void is not in good agreement with MD results [13]. (ii) The method implies Γ can be defined in the vicinity of the obstacle, but this is uncertain because the LT in the vicinity of the obstacle depends on the orientation of the line [19] and is not the same as that due to line curvature [20]. (iii) The force balance on a dislocation segment is not necessarily parallel to the segment direction [20]. A different method can be used to resolve these issues [21]. In the case of an obstacle in a periodic row, as shown in Figure 2b, imagine the obstacle to be equivalent to a dislocation segment of length w which experiences a stress preventing its motion. The obstacle experiences a force composed of two components: (a) the force τbw resulting from the applied stress τ, where w the shortest distance between dislocation arms touching the obstacle, as defined in Figure 2b, and (b) the force resulting from the curvature of the dislocation when the applied stress exceeds the friction stress, τ f , resulting from the lattice friction on the edge dislocation. The second component is given by lb(τ -τ f ), where l is the spacing between obstacles shown in Figure 2b. The following equation is imposed by the condition of force balance:

f eff w l w l τ τ τ -       + = 1 , ( 1
)
where τ eff is the 'effective' stress experienced by the obstacle when the applied stress is τ.

Mechanical equilibrium prevails as long as the obstacle counterbalances τ eff with an equal resistance in the opposite direction. At T = 0 K unpinning proceeds when τ eff exceeds the maximum obstacle resistance, o c τ , which can be deduced from Equation (1) when the applied stress reaches its maximum τ max . In most cases, l is much larger than w and this leads to the following approximation:

( )

f o c w l τ τ τ - = max .
(

It should be noted that Equation ( 2) is valid as long as the friction stress is independent of the dislocation character. Therefore, this method cannot be applied when a screw dislocation segment is formed during the interaction, which is the case for large void size [START_REF] Bacon | Dislocations in Solids[END_REF],

because of the large Peierls stress known to exist for the screw dislocation in iron at room temperature and below [8]. The idea of connecting unpinning to a critical effective stress rather than a critical force was already used successfully in dislocation dynamics simulations to predict strengthening by different mechanisms [22,23]. Here, we extend the method with the help of MD simulations to characterise void strength.

In the case of voids, the effective stress is associated with the stress felt by atoms on the void surface on both sides of the slip plane, tending to create the missing dislocation segment of length w that has been annihilated inside the void. The other term, o c τ , can be interpreted as the maximum resistance of the lattice to the creation of an edge dislocation at a curved free surface. o c τ should then be a characteristic material property and must be independent of the void size and simulation conditions. This point can be checked using Equation ( 2) if the critical shape of the dislocation is known for different void sizes and dislocation lengths.

An MS study for different void sizes was reported in [START_REF] Bacon | Dislocations in Solids[END_REF]. In order to check the validity of our method, we select the cases of relatively weak obstacles, i.e. D = 0.9, 1, 1. τ max obtained in [START_REF] Bacon | Dislocations in Solids[END_REF] with simulation boxes with length along the y axis equal to either 42 or 83 nm are plotted as circles or rectangles, respectively.

MD simulations results

Thermodynamic analysis at finite temperature is complex since thermal activation introduces a time-dependence into the interaction process and quasi-static treatment is no longer valid.

However, the first law of thermodynamics for energy conservation is still valid, but with the difference that the heat generated is not extracted from the system. This leads to an increase in temperature, ∆T, of the system.

The major contribution to this heat is dissipated work due to lattice friction [1]. It increases with decreasing temperature. The maximum value corresponds to the heat extracted from the system in MS simulation. The corresponding change in ∆T would be given by

        ∆ - = ∆ ∫ pot v U d V C T t γ γ τ 0 1 , (3)
where C v is a heat capacity at constant volume, dγ the increment of applied strain, γ t the total shear strain and ∆U pot the difference in the potential energy between the beginning and end of the simulation. The term in brackets in Equation ( 3) is the difference between the total mechanical energy provided to the system and the change in the potential energy, which corresponds exactly to the extracted heat. This heat amounts to 7.8 eV, which corresponds to ∆T < 1 K at the end of the simulation. Since this heat is extracted, the transformation is isothermal at 0 K. In MD simulations, the computed friction stress is so low (see observation (i) below) that it is difficult to assign a finite value to it. Even if the dissipated heat remains in the vicinity of the moving dislocation, the change in temperature in this region would be small. It is clear that this change can be neglected and the MD simulations, although irreversible, can be considered to be carried out at constant temperature and pressure, as a consequence of the relatively small friction stress on the edge dislocation.

To illustrate the effect of temperature, Figure 4a shows plots of stress versus strain for six different temperatures at the same applied strain rate and Figure 4b shows the four independent simulations at 300K and the same strain rate.

Fig. 4 : about here Several effects can be distinguished.

(i) τ f , which is 23 MPa at T = 0 K, decreases rapidly with increasing T, thereby confirming that the dissipated heat can be completely neglected.

(ii) τ max decreases with increasing T, which is the signature of thermal activation.

(iii) The pinning time decreases with increasing T.

(iv) The behaviour of the system is markedly stochastic, since the response is not identical in the independent simulations carried out at the same temperature.

Thus, unpinning occurs at different stress levels and different loading times. The latter is a specific feature of thermal activation. The amount of mechanical work necessary for unpinning is not the same in the different simulations and so a rigorous analysis should account for this. It will be seen in the following that the temperature effect and thermal activation can be analysed only on the mesoscopic scale. By 'mesoscopic' we mean the scale on which the atomic nature of the material is smoothed and replaced by other contiunuum models such as elasticity, elasto-plasticity, etc. For example, it is typically the scale used in dislocation dynamics simulations. At the mesoscopic scale, the atomic-level parameters should be replaced by homogenous variables. As the applied stress increases during straining from the reference strain γ = γ r , effective stress given by Equation ( 1) increases. This leads to relative displacement of atoms above and below the slip plane at the void surface, generating a local shear. A convenient way to describe the progress of the interaction is to introduce a scalar variable, ξ, representing the interaction coordinate, see Figure 5. However, it is difficult to assign a precise physical meaning to ξ in our case since atomic displacements are not restricted to the atomic rows belonging to the void surface.

Thermodynamic analyses

Free energy and the reaction coordinate

Fig. 5 : about here Following the discussion above, the transformation can be considered isothermal at constant pressure. The system potential is thus the Gibbs free energy, G, [24] and not the enthalpy. Figure 5 is a scheme of the interaction giving the evolution of the change ∆G(τ eff = 0) in G as a function of ξ in the absence of stress. Initially the system is in its ground state labelled I. To overcome the whole barrier, the free energy ∆G max should be provided in order to bring the system reversibly to its saddle state B. The derivative of the free Gibbs energy with respect to ξ gives the resistance stress τ lattice , i.e. ∂∆G(τ eff = 0)/∂ξ = bwτ lattice . A schematic representation of τ lattice is depicted in Figure 5. Under the effect of an applied effective stress τ eff , the system potential ∆G(τ eff ) evolves due to the work done by the effective stress, see Figure 5, and is now characterised by a new ground state labelled A at ξ

where τ lattice = -τ eff . Also, the saddle state becomes S and the energy barrier decreases to ∆G(τ eff ) because of the work done by τ eff during the transition from state A to state S. The system would stay at state A forever at 0 K. The unpinning of the dislocation occurs only when the effective stress is increased to its critical maximum value o c τ . At a temperature greater than 0 K, there is a non-zero probability for the missing free energy ∆G(τ eff ) to be supplied by the rest of the crystal, leading to the release of the dislocation.

It is important to note that the free energy profile and the corresponding reasoning are valid for quasi-static progress of the reaction. This, of course, is not the case during the thermal activation event, when the transition occurs in a time interval of the same order of magnitude as the period of atomic vibration [START_REF] Schoeck | Dislocation in Solids[END_REF]. The reaction coordinate has to be reconsidered correspondingly as a function of the inertial effects, which, to our knowledge, is 11 not covered by any theoretical or simulation investigation reported in the literature. This is the main problem of all approaches tending to estimate the activation energy directly from the Gibbs free energy profile [9,[START_REF] Seeger | [END_REF]27]. The approach proposed in this paper for determination of the activation energy allows this difficulty to be avoided.

Connection with statistical dynamics

Let us restrict our system to a small volume of atoms Ω surrounding the void. This volume is in thermal equilibrium with the rest of the crystal at temperature T. Each MD simulation can be considered as a sample path from an initial state I at time t = 0 to the saddle state S at t = ∆t, see Figure 5. An appropriate choice of these states facilitates the theoretical treatment. In our case, the natural choice for the initial and final states corresponds to the reference state γ r and the saddle state, respectively. The latter corresponds to the state of the maximum Gibbs free energy. Within every such path, work W is done by the effective stress on the volume Ω.

According to Jarzynsky [10,11] and Crooks [12], under the condition of microscopic reversibility, the increase of the system potential can be obtained by averaging the work done on the system over all sample paths:

( ) W G G G I S max β β - - = - = ∆ exp ln 1 , (4)
where β is the temperature parameter (kT) -1 and k is the Boltzmann constant. The brackets < > in Equation (4) denote the average over all paths. The average change in the system potential is, in general, less than the entire work done on the system, because of dissipation.

But Equation (4) stipulates that the modification of the system potential is a simple function of the average over all possible paths of the work done on the system. In our special case, all paths are characterised by the same initial and final states. The difference between G S and G I is therefore exactly equal to the height of the energy barrier ∆G max , see Figure 5.

Among the different paths in Equation ( 4), some of the statistical systems acquire the additional quantity of energy corresponding to the missing energy, ∆G = ∆G max -W, due to thermal fluctuations and therefore reach the saddle state. Once a system acquires ∆G from the surrounding crystal, it undergoes an irreversible transformation beyond the saddle state B.

This is why all paths that can be computed through MD simulations are restricted to only paths in which the system has acquired the activation energy. In contrast, the average in Equation ( 4) covers all possible paths and not only paths selected in MD simulations. This is the reason behind the difference between the expected change in energy ∆G ≡ W and the acquired value ∆G max .

Following the dynamical approach and equilibrium statistics [28], the probability of acquiring an additional energy ∆G is proportional to the Boltzmann factor exp(-β∆G). The exact jump rate ω [START_REF] Kocks | [END_REF] can then be given by ω = ν exp[-β ∆G], where ν is a frequency parameter to be discussed later. This rate equation stipulates that the probability of unpinning depends only on the instantaneous chance to receive ∆G from the rest of the crystal. During a time increment dt, the infinitesimal probability of acquiring ∆G is:

dp = ν exp (-β ∆G) dt (5)
The probability p, increasing with time, reflects the chance that every one of the statistical systems passes through the activated stress.

Activation energy under constant stress conditions

We investigate here the unpinning process in the simple condition of constant effective stress,

τ eff = τ c . Since τ c is lower than o c
τ , the system is in a stable equilibrium state called the ground state A in Figure 5. The system Ω is in thermal equilibrium with the rest of the crystal at finite temperature T. We assume that at t = 0 the volume Ω is in its ground state A, characterised by a given value of the thermodynamic potential. The dislocation unpins when thermal fluctuations provide the missing mechanical energy ∆G. Here we should note that (i) the energy barrier is assumed to be solely a function of stress and (ii) ∆G is a function of the effective stress defined in Equation (1), and not of the applied stress. This point is important since most treatments reported in the literature assume a dependency on the applied stress (for a review see [START_REF] Schoeck | Dislocation in Solids[END_REF]29]). Considering ∆G to have a dependency on the applied stress is not exact in the case of localised obstacles. It leads to unrealistic values for the activation volume, as will be seen later.

Since τ eff is constant, there is no reason for the Boltzmann factor not to be uniformly distributed over time. According to the transition state theory and dynamical theory [START_REF] Schoeck | Dislocation in Solids[END_REF]28], the average survival time, or the average time for the system to stay in the half space containing the ground state A, is simply the inverse, ω -1 , of the reaction rate. This survival time can be considered as the average incubation time of the unpinning process:

F
< ∆t > = ν -1 exp [β ∆G(τ c )]. (6)
Values of ∆t could be measured from MD simulations performed at constant applied stress and the Boltzmann factor, exp [-β ∆G(τ c )], then deduced from the average <∆t> using Equation ( 6). However, the simulations in the present work used constant applied strain rate loading and required the following treatment.

Definition of the critical stress in MD simulation performed at constant strain rate

For simulations performed under conditions of fixed applied strain rate, a definition has to be assigned to the critical stress τ c since the applied stress, and consequently the effective stress, varies with time. When τ eff is not constant, Equation ( 5) is only valid during an infinitesimal time increment dt. Here we have to distinguish between the maximum of the effective stress, τ max , related to the maximum of the applied stress through Equation ( 2), and the critical stress τ c , characteristic of the unpinning process. This issue was first addressed in [8], where a definition of τ c was provided. In MS simulations, the critical stress is equal to the maximum stress because the process is quasi-static, but at finite temperature the thermally-activated unpinning depends on time, which precludes the critical stress from being considered equal to the maximum stress. During loading from γ r in MD simulations, the effective stress increases with time, which leads to a rapid increase in the accumulated jump probability. The unpinning occurs after an incubation time ∆t. The critical stress is by definition the constant effective stress that would have provided the same jump probability during the same incubation time ∆t. This leads to the following expression: in which τ eff tends to τ max when t tends to ∆t. Since the Boltzmann factor is a monotonic increasing function of τ eff , it turns out from Equation ( 7) that τ c is always less than τ max . The difference between the two quantities depends on the temperature and the activation volume, which accounts for stress sensitivity of activation energy. In order to make progress in estimating the critical stress, we have to postulate a functional form of the activation energy.

[ ] ( ) ∫ ∆ ∆ - = ∆ - ∆ t eff c dt G G t 0 ) ( exp ) ( exp τ β τ β , (7)
As a first approximation, we can develop ∆G in the vicinity of τ max by ∆G = C -V * τ eff , where C is constant (not to be confused with the value of ∆G at zero stress) and V * the activation volume. It is important to underline here that our concept of the activation volume does not correspond to the apparent activation volume measured in experiment, V app . The connection with our definition would be V * = V app ×(w/l). Equation ( 7) allows us to deduce the critical stress:

t eff c V V ∆ = ) exp( ln 1 * * τ β β τ , ( 8
)
where the integral in Equation ( 7) is replaced by its average value over the incubation time.

The variation of V * in Equation ( 8) is relatively weak compared to the variation of τ eff . This allows us to deduce τ c using a rough estimate of V * . In our case, it can be set equal to wb 2 , which accounts for the interaction range of the process.

In order to illustrate the difference between τ max and τ c , consider the case where the effective stress increases linearly with time from 0 to τ max . The average in Equation ( 8) can then be derived analytically:

[ ] 1 ) exp( 1 ) exp( * * * - = max max eff V V V τ β β τ τ β . ( 9
)
It is easy to check that for typical values of β, V * and τ max , the exponential term in Equation ( 9) is dominant. The critical stress is then given by 10) confirms that τ c is less then τ max and shows clearly that the difference between them increases with temperature. On the other hand, when the temperature tends to zero, τ c tends to τ max . This asymptotic behaviour is expected since, in the absence of thermal activation, the unpinning proceeds only when the effective stress reaches the maximum stress.

kT V V kT c max * * max ln τ τ τ - = . (10)
Here we have to underline another difficulty in estimating the critical stress. In Equation ( 8), τ c provided by an MD simulation is only an estimate of the statistically representative critical stress. By repeating the MD simulation at the same strain rate and temperature, we obtain different estimates of τ c . However, averaging should be carried out in the probability space and not in the stress space because the stochastic feature of the activation has its origin in the probability of occupying a given state in the phase space. This probability is proportional to the Boltzmann factor, imposing an exponential dependency on the stress. The arithmetic (linear) average in the stress space is therefore not justified [9]. As a result, the value of τ c at a given strain rate and given temperature can be calculated from

n t eff c V V ∆ = ) exp( ln 1 * * τ β β τ , ( 11
)
where the inner average is to be carried out over time during every MD simulation, while the outer average is to be taken over the number n of repeated MD simulations. The validity of this approach is conditioned by the fact that the repeated MD simulations are independent.

The formula in Equation ( 11) is the same as that developed in [8]. The validity of the treatment in the latter work is thus established under an implicit condition: there is no correlation between the successive kink-pair nucleation events.

The treatment above has been used to calculate the effective stress in our simulations.

The results are shown in Figure 6, where it is seen how the maximum and the critical effective stress vary with temperature. However, in order to show the spread of values of τ c , we have used Equation ( 8) rather than Equation (11) and the different values obtained at each temperature are shown in the figure. It can be seen that the spread of values of τ c at a given temperature is quite small, so the convergence of the different estimates of τ c is sufficiently fast for our system. As expected, the difference between the stresses τ max and τ c tends to zero as T tends to 0 K. The difference increases rapidly with increasing T at low temperature, reaching an almost constant value of 0.7 GPa in the range 100 to 600 K. There are two reasons explaining this constant shift at high temperature. As can be seen in Eq. ( 11), the efficiency of the thermal activation is proportional to ln(τ max ). Also the slope of the function (kT)ln(kT)

decreases strongly with kT. It is thus easy to understand why the difference between τ max and τ c should saturate with increasing temperature. This difference, however, amounts to 20 to 25% of τ c . Consequently, it cannot be neglected and τ c should not be confused with τ max [9].

In other words, the unpinning process is not only determined by the maximum stress, but depends also on the time spent to reach this stress.

Calculation of the activation energy

The aim of this section is to provide a method for estimating the activation energy as a function of the critical stress. The unpinning process depends on three factors, namely the effective stress, the temperature and the deformation rate (through the incubation time), but only two are independent. Equation [START_REF] Marian | [END_REF] shows how these parameters are connected. Since in MD simulations the temperature is almost conserved and <∆t> can be estimated from different independent simulations, it is possible, in principle, to evaluate the activation energy using Equation ( 6), provided that the frequency parameter ν is known. This question has been addressed several times in the literature, for a review see [START_REF] Kocks | [END_REF][START_REF] Schoeck | Dislocation in Solids[END_REF]. The effect of the kinetic energy on the activation rate can be considered through the change in the vibrational free energy of the atomic row between state A and S [30]. The determination of ν for localised obstacles is complex and involves different characteristics [30]: strength, size, spread, etc. It also involves the vibrational spectrum of the bulk material. Large difference in the estimations of ν can be found in the investigations reported in the literature.

However, since the frequency parameter appears as a pre-exponential in Equation ( 6), only an order of magnitude estimate for its value is needed for the determination of the reaction rate. The commonly-accepted way of doing this is to scale the Debye frequency, ν D , by the length of the vibrating atomic row [31]. The frequency ν can be represented as the frequency of the transverse wave of length equal to that of the effective obstacle, which, in our case, is clearly equal to w, the smallest distance at which the dislocation is kept pinned, the vibration frequency of an atomic row of length w should be close to that of a dislocation segment of length w. By this procedure, we have ν ≈ ν D b/w. However, the vibrating atomic row belongs to the void surface, where the frequency spectrum is lower than that in the bulk [START_REF] Müller | Les surfaces solides : concepts et méthodes[END_REF], and the frequency should be corrected accordingly. As a first approximation, the slowing of atoms in the neighbourhood of the surface can be accounted for by assuming only half of the bulk frequency. For the situation where the stress is constant, Equation ( 6) can then be rewritten to estimate the activation energy:

      > ∆ < = ∆ w b t G D c 2 ln 1 ) ( ν β τ . ( 12
)
But for simulations at constant strain rate, the stress varies from τ f to τ max during loading. This is why the concept of the critical stress was needed in section 5.4 to provide the same jump probability as that measured during constant stress simulations. Equation (12) should thus be coupled to Equation (11). In other words, for every set of simulations performed at a given temperature, we have to evaluate independently the critical stress using Equation ( 11) and the activation energy using Equation (12). We then construct the two functions τ c (T) and ∆G(T) in order to deduce the desired function ∆G(τ c ).

The variation of ∆G with T obtained using Equation ( 12) is shown in Figure 7. The linear proportionality frequently observed in experiments on BCC and HCP metals such as iron [START_REF] Spitzig | [END_REF]34], zirconium [35] and titanium [36] is clearly confirmed in our results. This provides evidence for the validity of the two hypotheses that (i) the effect of stress and temperature on the attempt frequency is negligible and (ii) the profile of the energy barrier, i.e. the mechanical work necessary to overcome the barrier without the help of thermal activation, is independent of temperature. The linearity can be described by the relation ∆G =

CkT, where C is a constant close to 8.1. The values C = 8.1 and 6.6 correspond to & γ = 10 6 and 5x10 6 s -1 , respectively. The change in activation energy confirms our assumptions. Increasing the strain rate decreases ∆G because a faster increase of the applied stress leads to a decrease in the incubation time.

Combining data from Figure 6 for τ c with that for ∆G(τ c ) from Figure 7 enables us to determine the dependence of the activation energy on the critical stress, as shown in Figure 8.

Over most of the range, the dependence is almost linear and consistent with a constant value of the activation volume via the relation ∆G(τ c ) ∝ -V * τ c , which was used to deduce Equation ( 6). The data give V * = 3.6 b 3 , which is close to our estimate made for the evaluation the critical stress. This ensures that the procedure for the estimation of the critical stress is compatible with that used for the evaluation of the activation energy. Beyond the linear region of the plot in Figure 8, ∆G(τ c ) tends to zero when τ c approaches the value found for T = 0 K. Fig. 8 : about here

The values of the activation energy presented above were measured for critical stresses larger than 2.5 GPa. The activation volume does not appear to increase as rapidly as expected when the critical stress falls below this value and tends to zero. Creation of the missing dislocation segment at the void surface should be accompanied, in principle, by a long-range perturbation of the lattice. Although the corresponding mechanical work would still be finite, the activation energy should increase rapidly for low stress [START_REF] Kocks | [END_REF]. In order to investigate the low stress regime, MD simulations should be carried out at higher temperature and/or much lower strain rates. Unfortunately, this is not feasible with current computer resources.

Conclusions

It has been shown that a criterion for the obstacle strength of voids, expressed in terms of a critical local effective stress, can be derived from molecular statics calculations carried out to simulate a crystal at temperature T = 0 K. The critical effective stress was found to be independent of the void size, the dislocation length and other simulation parameters. The efficiency of thermal activation in decreasing the critical effective stress is found to be a They represent the change of system energy at zero effective stress, the system energy at a given effective stress and the effective stress needed to push the interaction to a given ξ. We simplified the plots by considering the friction stress to be zero.

  interaction strength in the case of finite size obstacles. The classical picture

1 .

 1 Dependence of stress and energy on applied strain

5 and 2

 52 nm, for which unpinning is controlled by dislocation curvature. Data for larger voids cannot be used because MS simulations show that straight screw dipoles are generated at the obstacle before unpinning, which connects the unpinning process to the mobility of screw segments. For the different sizes considered, the values of w are taken to be equal to 4a 3 a is the lattice parameter. These values correspond to the integer number of interplanar spacings along ] 2 11 [ that are closest to the observed critical width. According to Equation2, a plot of τ max versus w/l should give a straight line with intercept τ f (approx. 23 MPa) on the τ max axis. This is tested in Figure3, where MS results for

Fig. 3 :

 3 Fig. 3 : about here

Fig. 6

 6 Fig. 6 : about here

17 see Figure 2 .

 172 Since the result of activation is formation of the missing dislocation segment,

Fig. 7 :

 7 Fig. 7 : about here

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2. (a) Conventional picture of mechanical equilibrium for a point obstacle in the line tension approximation. (b) Equilibrium around a finite size obstacle.

Figure 3 .

 3 Figure 3. The maximum applied stress as a function the ratio w/l for voids (see text for more detail). The dashed line is the best linear fit.

Figure 4 .

 4 Figure 4. Stress-strain plots for (a) six temperatures and (b) four simulations for T = 300 K. & γ = 10 6 s -1 in all cases.

Figure 5 .

 5 Figure 5. Schematic plots of interaction profiles as a function of the interaction coordinate, ξ.

Figure 6 .

 6 Figure 6. Dependence of the maximum stress, τ max , and the critical effective stress, τ c , on temperature for all MD simulations. The values of τ c were calculated from Equation (9).

Figure 7 .

 7 Figure 7. Variation of ∆G with T obtained using Equation (13).

Figure 8 .

 8 Figure 8. Variation of the activation energy as a function of the critical stress, τ c , obtained from Figures 6 and 7.

  Fig.5

Fig. 7 Page

 7 Fig.7

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgments

This work was partially supported by the Division of Materials Sciences and Engineering and the Office of Basic Energy Sciences, U.S. Department of Energy, under contract with UT-Battelle, LLC. It was also supported by the European project PERFECT ( FI60-CT-2003-208840 and by grant GR/S81162/01 from the UK Engineering and Physical Sciences Research Council.