Alexander V Barashev
email: a.barashev@liv.ac.uk

Stanislav I Golubov

Yuri N Osetsky

R E Stoller

Dissociation of migrating particle from trap with long-range interaction field

Keywords: diffusion, kinetics of reactions user supplied): reaction-diffusion theory, defect complexes, dissociation lattice defects, reaction-diffusion theory, defect complexes, dissociation

Dissociation of migrating particle from trap with long-range interaction field

Introduction

Solving various problems of diffusion-reaction kinetics requires knowledge of the time a migrating particle spends trapped by lattice defects before dissociation. This information allows considering longer time and length scales using a rate theory [1] or a kinetic Monte Carlo approach, which treats complexes as single entities. These dissociation times have been or can readily be obtained only in simple cases, e.g. for vacancy-solute pairs with shortrange bonding. In more complicated cases, the solution is generally unknown. A particular example of such a case is that of a vacancy executing three-dimensional (3-D) random walk, which is trapped inside a Cu precipitate in an iron matrix (see [2][3][4][5] for molecular dynamics results). Another example is that of a one-dimensionally migrating cluster of self-interstitial atoms trapped in a long-range field of an edge dislocation (see, e.g. [6] for observations in Ni and Mo). Yet another example is that of an interstitial cluster trapped between Cu precipitates in Fe; it can penetrate inside one of them and this requires overcoming an energy barrier (see [4,5] for molecular dynamics results). The importance of the processes mentioned above in microstructure evolution and changes of mechanical properties, especially under irradiation conditions, is discussed in [1][2][3][4][5] for Cu precipitation and, in [6,7],

for dislocation decoration with interstitial clusters.

In this paper, we first derive an equation for the mean dissociation time for a well of arbitrary shape in section 2. We start with a simpler case of one-dimensionally diffusing particles in section 2.1, and then generalise to higher-D diffusion in section 2.2. Some particular potential energy wells and limiting cases are considered in section 3. Monte Carlo calculations of the probability distribution of dissociation times are presented in section 4 for two particular 1-D wells. In section 5, we generalise the problem to the case of particle exchange between spatially separated wells for two particular cases: vacancy evaporation from void and interaction of interstitial clusters with dislocation. The conclusions are drawn in section 6.

Mean dissociation time from well

One-dimensional diffusion

Consider a particle that hops on a lattice containing a potential well U(x i ) (interaction energy in equilibrium positions i, see Figure 1a) in the diffusion direction and 'free' states. The 'free' states are those, where U (x i ) = 0 and all the migration barriers are the same as in perfect crystal. The mean time-delay between jumps in these states is denoted by τ . A dissociation event occurs when particle jumps from the well to a 'free' state and an association event is a reverse jump. The dissociation time τ diss is the mean time from association to dissociation. (This definition is discussed further in sections 4. 3 and 5) [Insert Figure 1 Consider a well and one adjacent 'free' state in periodic boundary conditions. In this system, a dissociation event is followed by an association event, hence for each time that particle spends time τ diss on average in the well region, it spends time τ in the 'free' state.

Hence, the probability to find particle in the 'free' state is equal to

p free = τ τ + τ diss . ( 1
)
The same probability can be written through the partition function of the canonical ensemble

Z = 1 + e -βU ( x i ) well ∑
, where β = 1 / k B T , k B is the Boltzmann constant, T is the temperature and unity corresponds to the 'free' state, where U = 0 , (see, e.g. [8]) as:

p free = 1+ e -βU ( x i ) i=1 Ν ∑       -1 . ( 2
)
In this equation and further in the text, the summation is taken over all Ν equilibrium states of the well. From equations ( 1) and ( 2) one obtains

τ diss = τ e -βU ( x i ) i=1 Ν ∑ = τΝ e -βU ( x i ) well , (3)
where the brackets denote averaging over sites representing equilibrium states of the well.

Higher-dimensional diffusion

The analysis can readily be extended to two and three dimensionally diffusing particles.

Consider a well U(r i ) of Ν states surrounded by a shell of adjacent 'free' states, m in total, each separated from the well by one jump distance. Assume that a particle in a 'free' state jumps successfully only to the well. This assumption does not influence the results but makes the derivation easier. It affects, however, the times spent by a particle in 'free' states. If the mean value of these times, τ free , is known, the probability of finding the particle in a 'free'

τ free τ free + τ diss = m m + e -βU (r i ) i=1 Ν ∑ . (4)
Hence, the dissociation time is given by

τ diss = τ free m Ν e -βU (r i ) well . ( 5
)
Now let us find τ free . Under the assumption that a particle in a 'free' state jumps successfully only towards the well, the time spent on average by the particle in a 'free' state i before jumping back to the well is higher than τ by the ratio of the total number of jump directions z to the number of directions towards the well, n i : τ free i = τ z / n i . The mean time is the sum of these times weighted with relative frequency of visiting these states, p free i , which is equal to the ratio of the number of directions from the well to this particular site, n i , to the total number of directions from the well to 'free' sites,

n i i =1 m ∑ . Hence p free i = n i / n i i =1 m ∑ and τ free = τ free i p free i i=1 m ∑ = τ z / n free , (6)
where the brackets denote averaging over adjacent 'free' states. Finally, by substituting equation (6) into equation ( 5), one obtains

τ diss = zΝ m n free τ e -βU (r i ) well . (7)
This is a general solution valid for any dimensionality of diffusion. In the case of 1-D diffusion, when the total number of jump directions is z = 2 , the number of directions from a We emphasise that equation (7) depends on the energy of equilibrium states (the exponent is averaged over these states) but is independent of the saddle point energy profile within the well. Note that τ depends on the energy barrier in the perfect part of the crystal only. (Obviously, the probability distribution of dissociation times does depend on the saddle point energy and this is considered below in section 4.) Equation ( 7) is exact within the assumptions used, i.e. it gives the mean time a diffusing particle spends within a well region between association and dissociation events, which are defined as particle crossings of a dividing surface between the potential well and the rest of the system. It should be emphasised that equation (7) has the meaning of the time when the number of particle-potential well complexes reduces by e times only for exponential distribution of dissociation times. This is thus a necessary condition for the use of equation (7) together with the rate equations without memory effects. This in turn requires separation of the time scale of dissociation from other fast scales, e.g. time for crossing the barrier and equilibration time of a particle within the well (see e.g. [START_REF] Hanggi | [END_REF], p.257). In particular, to ensure separation of the time scales, the potential well must be deep enough compared with the thermal energy, β -1 .

3. Mean dissociation times for particular wells

1-D square well

For a square well and coordinate x (in units of jump length a )

U (x) = -E, x min ≤ x ≤ x max , 0, x < x min , x > x max .      (8)
Equation (3) yields where Ν = x maxx min is the total number of equilibrium positions inside the well. Hence the frequency of dissociation events is given by

τ diss = Ντ e β E , (9)
ν diss = τ diss -1 = 2 Ν       e -β E 2τ . ( 10
)
The physical significance of equation ( 10) can be understood as follows. The bracketed term stands for the probability to find a particle at one of the two well edges, τ -1 is the total jump frequency in perfect crystal in both directions, ½ describes the probability to jump towards the barrier and the exponential factor accounts for the decrease of the probability of successful jump due to the barrier.

1-D triangular well

A particular case of a triangular shape well:

U (x) = -E(1-x / x min ), x min ≤ x < 0, -E(1-x / x max ), 0 ≤ x ≤ x max      (11)
is a more realistic approximation to real situations with no single energy. In this case:

τ tri = 2τ e β E 1-e -β E (1+ x max -1 )     1-e -β E / x max . ( 12
)
It can readily be shown that, for deep enough triangular wells, such that e βE >> 1 , the dissociation time is smaller than for a square well of the same depth and width by the ratio of the thermal energy to the well depth. This is due to a smaller effective binding energy.

Void

A void can be considered as a well for vacancies. Since vacancies are always at the void surface, Ν = m in equation (7). Hence, for a large void of radius r 0 , such that n free = z / 2 , successful jumps are only in the direction away from the void. The evaporation frequency is, hence, ν diss = e -β E / 2τ . In equilibrium conditions this frequency equals to that of the void- vacancy association events: C V e / 2τ , hence the equilibrium vacancy concentration near the void is C V e = e -β E (see section 5 for derivation). This equilibrium concentration and equation ( 7) can be used to calculate the evaporation rate using the diffusion-reaction theory (see section 5 below).

Vacancy in solute precipitate

For a vacancy and a large solute precipitate of radius r 0 with a constant interaction energy between the vacancy and any site within the precipitate U = -E : n free = z / 2 , Ν = 4πr 0 3 / 3Ω , where Ω is the atomic volume, m = 4πr 0 2 a / 2Ω , where the surface layer is assumed to be of thickness a / 2 , and, hence, ν diss ≈ e -β E a / r 0 τ . The ratio a / r 0 accounts for the probability to find a vacancy near the surface. Hence, the frequency of dissociation events is equal to the vacancy concentration near the surface times the frequency of jumps away from the precipitate, which seems to be an obvious description in such a simple case.

To take into account the interaction energy profile U(x i ) within the precipitate, one should use e -βU (r i ) well -1

for the probability of successful jump in this case. This is useful in calculations of an enhanced vacancy concentration inside precipitate.

For example, in conditions typical for thermal ageing, when the concentration in the matrix outside a precipitate is kept constant and equal to the equilibrium vacancy concentration, C V e , at given temperature, the enhanced mean vacancy concentration in the precipitate C in can be obtained using the detailed-balance condition at the precipitate surface. The frequency with which vacancies enter the precipitate is C V e m n free / zτ , while the frequency with which vacancies leave the precipitate is C in Ν / τ diss . The equality of these frequencies implies that 3) and ( 7), is given by the probability density function of dissociation times, dP(t) / dt . Here, P(t) is the probability that a dissociation event occurs before time t . For a random process without memory, such as the radioactive decay:

1 -P(t) = exp(-t / τ diss ) (14)
and the probability density function is given by

dP(t) / dt = exp(-t / τ diss ) / τ diss . (15)
It will be shown below that this is valid in the cases characterised by a unique binding energy such as a single-site trap, but not for spatially distributed wells.

To study the influence of the well shape on the distribution of dissociation times, we performed Monte Carlo calculations of particle diffusion in the same system as used for derivation of the mean dissociation times in section 2, i.e. it consisted of a well and one 'free' state with periodic boundary conditions. Two particular cases are considered below, namely 1-D square and triangular wells (see sections 3.1 and 3.2). In the case of a square well, the jump frequencies of the diffusing particle in forward and backward directions, ν + and ν -, were equal to those in perfect crystal:

ν + (x) = ν -(x) = 2τ
( ) -1 for all but two frequencies for two border sites, where ν right border

+ = ν left border - = e -β E / τ .
In the case of triangular wells, the frequencies were calculated as ν ± (x) = e m β∆U ± / τ , where

∆U ± = ± U(x ± 1) -U(x) [ ]/ 2 are
the changes of the migration barriers due to the interaction energy (see Figure 2). The probabilities that particle jumps forward or backward are defined by ratios ν + / ν total and ν -/ ν total , where ν total = ν + + ν -is the total jump frequency. The value < t >= ν total -1 is the mean waiting time before a jump, which is generally site dependent. The value τ normally obeys the Arrhenius relationship for the temperature dependence with activation energy equal to the saddle point energy in perfect crystal ( E 0 in Figure 2). In most of our calculations, τ was a constant chosen to be the unit of time. Some calculations revealing influence of the saddle-point energy profile on the distribution function, thus with τ dependent on the position inside the well, are presented in section 4.4. The calculations were performed until a specified number of dissociation events occurred. The distribution of dissociation times was obtained as the number of dissociations within a specific time window (a bin of the distribution) divided by the total number of events and the bin width. In our calculations, the bin width was taken as 0.1τ diss and calculated using equations ( 9) and ( 12) for the square and triangular wells, respectively. The main purpose of calculations was to compare the probabilities calculated with those described by equation (15). We use a dimensionless variable t / τ diss and function τ diss dP(t) / dt to represent the results, since they are independent of the mean dissociation time.

[Insert Figure 2 about here]

Results for 1-D square wells

The MC calculations were performed for wells of different width (50, 100 and 200 one jump distances) and depth ( βE =3.87 and 7.74, which correspond to E =0.1 and 0.2 eV, respectively, at T =300K) and ~500,000 dissociation events were accumulated in each case.

The range of well widths was chosen with a particular interest in the migration of a selfinterstitial atom cluster in the long-range stress-field of an edge dislocation [6,7]. The calculations reproduced equation [START_REF] Hanggi | [END_REF] with high accuracy, which was one of the tests for the code. The main aim of the following calculations was to study the deviation of dP(t) / dt from equation (15). The results are presented in Figure 3. Typical error bars are shown for one calculation (the upper parts of errors bars are omitted in order to keep presentation clear).

As can be seen, for wells deep and/or narrow enough (e.g. for βE =7.74 or Ν =50), the probability density approaches equation (15), which is characteristic of a random process where the mean dissociation time is not significantly longer than τ diff = Ν 2 τ / 2 , which is the time corresponding to the diffusion length equal to the half of the well width l = Νa / 2 , where a is the jump distance. The diffusion length, i.e. the mean-free path, is defined from equation l 2 = 2Dτ diff , where D = a 2 / 2τ is the diffusion coefficient in 1-D. The last column of the legend in Figure 3 contains necessary ratios τ diss / τ diff . We note that wells for which τ diss ≈ τ diff can hardly be treated as traps at all. Thus, we conclude that for square wells strong enough to be considered as traps, i.e. with τ diss >> τ diff , the distribution under interest is described reasonably well by the exponential dependence. Note that the condition τ diss >> τ diff is one of the basic ones required for separation of the time scales (see section 2).

[Insert Figure 3 about here]

Results for 1-D triangular wells

Calculations similar to those described above for a square well were performed for triangular symmetrical wells, i.e. when x min = x max in equation (11). The statistics were accumulated for about 1,000,000 events. Equation (12) for the mean dissociation time was reproduced within reasonable accuracy, usually about 1%. The distribution function of dissociation times is presented in Figure 4. As can be seen, the distribution deviates from the exponential dependence much more significantly than for square wells. Similar to square wells, the deviation is stronger for shorter times and smaller ratio τ diss / τ diff . This is because of a high fraction of dissociation events occurring before particle visited the deep region of the well. It is evident that the higher the ratio τ diss / τ diff , the lesser the contribution of shallow regions to the mean time. Unlike that for the square well, the probability of short time dissociation in a triangular well is always high, even for deep and narrow wells, and the distribution at short times has quite different slope in logarithmic scale. This feature is due to geometry of the well, where shorter times correspond to smaller energy. At long times the distribution becomes exponential but with different slope, as shown in the comparison with exponential function plotted in Figure 4. Figure 5 shows the distribution of 1 -P(t) , i.e. the probability that a dissociation event does not occur before time t, for triangular wells of different width and depth. The error bars were very small for calculations of this type, not visible on the graph. One calculation for a square well (labelled 'Square') is also shown and it demonstrates a dependence which is close to exponential decay. In contrast, the distributions for triangular wells deviate from simple exponential dependence. The time dependence exhibits different behaviour at short and long times. The function at short times is steep and decreases quickly, while an exponential dependence evolves at long times. The curve for βE =3.87 and Ν =100, which is labelled 'x=0', is the distribution of escape times from the bottom of the triangular well. That is, the trajectories of all the diffusing particles in this calculation started from x = 0 corresponding to the deepest position in the well. In this calculation the mean dissociation time was about 15200τ and more than an order of magnitude longer than that given by equation ( 12) and corresponding calculations, which is 1210τ . We also note that this time was about twice higher than the corresponding value of 7230τ for the square well of the same depth and width when all trajectories started from 'x=0', i.e. the middle of the well. This does not seem trivial, since the simple average depth of the former well is smaller.

As seen from Figure 5, this calculation also shows the exponential decay. This indicates that the deviations from an exponential relationship are due to trajectories, which do not reach the bottom of the well. For such trajectories the dissociation times are short, so that the time scales are not well separated (see end of section 2). In addition, we make an important observation. When the time is equal to the mean dissociation time, the probability for the triangular wells is not equal to e -1 as expected for a random process, such as radioactive decay, described by equation ( 14). This means that the dissociation events are not randomly distributed, i.e. their time distribution does not obey equation ( 14), and hence actual distribution function is required to characterise the process.

It seems important to mention that, strictly speaking, such a process cannot be treated by chemical reaction equations, where the dissociation rate is the reciprocal of the mean dissociation time. These equations predict exponential decay for the trap concentration, thus describe random process. In reference [10] we show that by redefining the dividing surface

(

Influence of saddle-point energy profile on time distribution

As has already been mentioned in section 2.2, the mean dissociation time is independent of the saddle point energy profile within the well, and this follows from the general equation (7).

We verified this statement by calculations for a triangular potential well where all saddlepoint energies were taken to be the same and equal to those in perfect crystal. An example of such a well is shown in Figure 1b. In the calculations performed, the time delay was defined according to the actual saddle-point energy by including an additional Arrhenius factor. The results are presented in Figure 4 in the curve labelled 'SP' for βE =3.87 and Ν =100. The mean dissociation time was the same as in corresponding calculation above for triangular well (see section 4.3) and given by equation ( 12) within accuracy of less than 2%. The probability distribution was different, however, with the probabilities of shorter-time dissociations increased. This is because there is no longer preference for a diffusing particle to jump towards deeper well regions. In other words, in this calculation, the well represents a region where only jump frequencies, but not the probabilities of different jumps, are affected.

Similar effect was studied in [11] for diffusion of vacancy and Cu atoms in a dilute

Fe-Cu alloy. The analysis showed that the diffusion mechanism of Cu atoms, namely the crossover from 'exchange' (when a vacancy executes only exchange jumps with Cu atom) to 'drag' (when the vacancy co-migrates with Cu atom) mechanisms depends critically on the saddle-point energy profile of the vacancy around Cu atom (see [11], especially section 3.2).

Exchange frequency between wells

In the preceding sections we used the term 'dissociation', which can be inappropriate in some cases where the term 'escape' would be a better word. An important example is a dilute solution of traps. In this case, a migrating defect escaping from a trap has a high probability of returning to the same trap after a short-time migration in pure crystal. In diffusion-reaction kinetics, however, we are concerned mainly with events that result in complete dissociation in the sense that the memory is lost and the next reaction will occur with high probability at In other words, the quantity of interest is the frequency with which diffusing objects are exchanged between wells. This is the origin for the Eyring transmission factor in the reaction rate theory [START_REF] Hanggi | [END_REF]. A clear distinction between escape and dissociation events can be seen in the example of vacancy evaporation from voids. The evaporation takes place at the void surface; hence, the total escape frequency is proportional to the surface area, i.e. the void radius squared. The well-known result of the diffusion-reaction theory gives proportionality to the void radius. Obviously and as demonstrated below, this difference is due to correlated absorption of vacancies by the same void they are emitted from. This indicates that accounting for correlated events can change even qualitative description of the process. Note also that the events contributed to the short-time non-exponential part in the distribution function of dissociation times, see section 4.3, can also be viewed as correlated events and thus part of the same problem. Below we consider this problem for two specific cases:

vacancy evaporation from a void and de-trapping of interstitial cluster from dislocation.

Vacancy evaporation from void

Consider a void of radius r 0 , which emits τ diss -1 vacancies per second per surface site, in a spherical coordinate system. Vacancies migrate three-dimensionally with the diffusion coefficient D = a 2 / 6τ . The diffusion equation for vacancy concentration C is

∇ 2 C = 0 . (16)
To calculate the number of vacancies emitted from the void that reach some distance R from the void surface we use absorbing boundary conditions at this distance

C(R) = 0 . (17)
One more boundary condition must specify the vacancy-void interaction process.

Assuming that vacancies are absorbed by the void, which is a realistic scenario, the vacancy concentration at one jump distance a from the surface can be written as The left-hand side of the equation describes the frequency with which vacancies leave the site. The first term on the right-hand side stands for the production of vacancies due to evaporation from the void. The last term on the right-hand side accounts for vacancies coming to this site from the sites further way from the void surface. After representing the latter term using Taylor series, in the limit of r 0 >> a , the boundary condition, equation (18), assumes the following form

C(r 0 + a) / τ = ν diss + C(r 0 + 2a) / 2τ . (18)
C(r 0 ) = 2τν diss + ∇C(r 0 )a . ( 19
)
Using this condition, one finds the concentration to be equal to

C = 2τν diss r -1 -R -1 r 0 -1 -R -1 . ( 20
)
It can readily be estimated using the last two equations that the gradient of concentration in equation ( 19) is smaller than the other terms by a factor of a / r 0 and does not contribute to equation ( 20). This means that most vacancies emitted from the void return back to it. As a result, the equilibrium condition for the concentration near the void surface is defined by the equality of frequencies of evaporation and jumps back to the surface and is not affected by the flux of vacancies away from the surface. The vacancy equilibrium concentration at the void surface is readily obtained from equation (20) as

C V e = C(r 0 ) = 2τν diss .
The total number of vacancies passing through a spherical surface of radius R and area S = 4π R 2 per unit time is equal to

R V a = - SD Ω ∇C(R) = DC V e Ω 4π r 0 1 -r 0 / R , (21)
where superscript 'a' denotes absorbing boundary. There are three points to be made. First, equation (21) becomes independent of the distance R from the surface, when R >> r 0 . Thus, vacancies reaching this distance lose their memory and can be counted as dissociated from the void. Second, despite the fact that the total vacancy emission frequency is proportional to the void surface area, the total vacancy flux far away from surface is proportional to the void radius. This is a well-known result of the reaction-diffusion theory (see, e.g. [12]). Third, as can be seen from equation ( 21), significant deviation from the proportionality to the void radius occurs at distances of the order of void radius.

As discussed above most vacancies emitted return to the void. The fraction of vacancies which do not return is equal to the ratio of the frequency defined by equation ( 21) and ~4π r 0 2 ν diss / a 2 , which is for the total frequency of vacancy emission. It is thus equal to a / r 0 . The same result can be demonstrated considering another, although unrealistic, scenario in which vacancies are reflected by the voids. In this case all vacancies produced contribute to the flux far away from the void. Indeed, in this case, we should multiply the left hand side of equation ( 18) by ½ to account for the fraction of unsuccessful jumps towards the void surface. Then, only the gradient and source terms remain in the boundary condition, equation

:

∇C(r 0 )a = -2τν diss ( 22
)
and the vacancy concentration becomes

C = 2τν diss r 0 2 a 1 r - 1 R       . (23)
The quantity defined by equation ( 21) modified for the case of reflecting boundary conditions (superscript 'r') is now proportional to the surface area:

R V r = 4πr 0 2 Ω 2Dτν diss a = r 0 a       R V a , (24)
and is higher than that for absorbing conditions by a factor of r 0 / a .

We also note that the first non-vanishing correction to the proportionality of the vacancy flux to the void radius is positive and proportional to the void radius squared, see equation ( 21), where r 0 1 -r 0 / R ( )

-1 ≈ r 0 + r 0 2 / R . This is the same result as obtained previously e.g. by Göselle [13] considering void capture efficiency. Thus, with increasing volume fraction more and more vacancies become absorbed at other voids and the shows the right tendency. Similar transition of the radius dependence of the reaction rate from R to R 2 can be found in [14] for precipitation kinetics for different absorption abilities of the solute precipitate-matrix interface.

Interstitial cluster exchange between dislocations

Consider a random spatial arrangement of dislocations and interstitial clusters migrating onedimensionally with the diffusion coefficient D = a 2 / 2τ . This problem is similar to that considered above in section 5.1. It is readily obtained that, with absorbing boundary conditions at the well boundary at x = 0 , equation ( 18), and at distance L from the well, the cluster concentration is described by the following equation

C = 2τν diss 1 - x L       . ( 25
)
The distance L here can be interpreted as the mean distance between dislocations in 1-D and will be discussed further below. In this case, other dislocations absorb clusters passing this distance and, hence, the memory is lost. This process can be considered as cluster exchange between wells represented by dislocations. The total flux of clusters emitted from the well at distance L is

J a = -D∇C(L) = D Ω 2τν diss L . ( 26
)
We note that, in contrast to 3-D diffusion, the total flux for 1-D diffusion always depends on the distance; namely, it decreases with increasing distance (decreasing dislocation density).

Similar to the problem of vacancy evaporation from a void considered above, one can show that, for reflecting boundary conditions, equation ( 22), the vacancy concentration and the flux of clusters are both L / a times higher. Hence, the fraction of clusters, which do not return to the well after escaping from the potential well, and, hence, the frequency with which clusters are exchanged between dislocations, f is equal to a / L :

Now let us make some estimates of the effect. The mean distance between randomly arranged dislocations in 1-D is described by the dislocation density ρ D and capture radius d as L = 2 / πρ D d (see, e.g. [15]). The jump distance in bcc lattice is the distance between atoms in <111> direction: a = 3a 0 / 2 , where a 0 is the lattice parameter. Hence, the frequency of cluster exchange between dislocations is

f = 3πν diss ρ D da 0 / 4 . (29)
For typical dislocation density ρ D =10 12 m -2 , the capture radius d =100 nm estimated in [6] and a 0 =0.3 nm, the frequency is smaller than the rate of cluster emission from the well by about five orders of magnitude: f ≈ 4 × 10 -5 ν diss .

Thus, the exchange frequency of the diffusing particles between spatially separated wells may be many orders of magnitude smaller than for frequency of dissociations from the well, and this is due to correlated recapture of the particles by the same well. And the effect depends strongly on the trap density.

Conclusions

1. An equation for the mean dissociation time of a migrating particle from a trap has been derived. It is independent of the saddle point energy profile within the well.

2. Generally, the distribution of dissociation times deviates from exponential function, especially in the region of small dissociation times. In such cases the probability function at the mean time may differ significantly from 1e -1 , which is for a random process. The effect depends on the well shape (is stronger for shallower and/or wider wells) and on the saddle point energy profile within the well.

3. The exchange frequency of the diffusing particles between spatially separated wells is generally many orders of magnitude smaller than the frequency of dissociations from the well, and this is due to correlated recapture of diffusing particles by the same well.

  written in two ways: as the ratio of times (left-hand side of the following equation) and via the partition function (right-hand side):

5 '

 5 free' state to the well is n = 1 and the number of 'free' states adjacent to the well is m = 2 , equation (7) reduces to equation (3).

  binding energy U E -= , one obtains τ diss ≈ 2τ e β E . The factor of two in this equation accounts for an increase of the time spent by a vacancy in the void, because

1 .

 1 Definitions and Monte Carlo schemeA more detailed description of the stochastic process of dissociation than by the statistical average value of dissociation time, equations (

  For shallow and/or wide wells, the distribution deviates from exponential dependence, especially at small times. More specifically, the deviation occurs for such wells,

11 [

 11 Insert Figure4about here]

12 [

 12 boundary of the well) the fast time scales can be eliminated and the exponential distribution of the dissociation times thus restored, Insert Figure5about here]

  void radius squared would be restored. The first correction term just

  diss a / L .

Figure. 1 .

 1 Figure captions

Figure. 2 .

 2 Figure. 2. Schematic diagram illustrating calculation of energy barriers.

Figure. 3 .

 3 Figure. 3. Probability distribution of dissociation times of a migrating particle from wells of different depth E and width Ν . The typical error bars are shown for one set of calculations (the upper parts of errors bars are omitted in order to keep presentation clear).

Figure. 4 .

 4 Figure.4. Same as in Figure1but for triangular wells. One curve, labelled 'SP' has been calculated for the energy profile where all saddle-point energies are equal.

Figure. 5 . 1 Figure. 1 . 2 Figure. 2 . 3 Figure. 3 .

 5112233 Figure.5. Probability that a dissociation event does not occur before time t for the wells of triangular shape of different width and depth. One curve labelled as 'Square' was calculated for square well. The curve labelled as 'x=0' shows the distribution of escape times from the bottom of the triangular well.

4 Figure. 4 . 5 Figure. 5 .

 4455 Figure. 4. Same as in Figure 1 but for triangular wells. One curve, labelled 'SP' has been calculated for the energy profile where all saddle-point energies are equal.

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

Acknowledgments A.V.B. acknowledges a research grant from the UK Engineering and Physical Sciences Research Council. Research at ORNL was sponsored by the Division on Materials Sciences and Engineering (R.E.S. and Y.N.O.) and the Office of Fusion Energy Sciences (S.I.G.), U. S. Department of Energy, under contract no. DE-AC05-00OR22725 with UT-Battelle, LLC.


Angus J Wilkinson
email: angus.wilkinson@materials.ox.ac.uk

David Randman
email: d.randman@sheffield.ac.uk

Determination of elastic strain fields and geometrically necessary dislocation distributions near nanoindents using electron back scatter diffraction

Keywords: dislocation mechanics, dislocations, electron diffraction, electron microscopy, mechanics of materials, nanoindentation

Determination of elastic strain fields and geometrically necessary dislocation distributions near nanoindents using electron back scatter diffraction

Angus Wilkinson, David

Introduction

Indentation has been used for a long time as a convenient means of probing the mechanical properties of materials. The introduction of instrumented indenters has significantly increased the flexibility of indentation testing and interest in testing thin films and small material volumes has generated a large user base for nanoindenters. Indentation testing has been applied to a wide range of materials systems to generate mechanical property data including hardness, elastic modulus, film delamination, residual stresses, and fracture toughness. One of the fundamental aspects of indentation testing is the well known indentation size effect (ISE). Although hardness values remain constant for large indents made with conical or pyramidal indenters, as the indentation depth is decreased below a micron or so there is a marked increase in the measured hardness. This ISE has been demonstrated in many materials including for example silver [1], tungsten and molybdenum [2], iridium [3] and is not an artefact of sample preparation. Nix and Gao [4] have given a convincing explanation of the effect in terms of geometrically necessary dislocations (GNDs) that are generated as the indent is formed. Some modifications have been made to the original analysis to account for the finite volume around the indent in which the GNDs are accommodated. However, in each case the GNDs are accommodated in smaller volumes of material for shallow indents leading to an increase GND density. Nix and Gao [4] showed that the average GND density near the indenter should be proportional to the reciprocal of the indenter depth. At sufficiently shallow indents the GND density can dominate over the statistically stored dislocations and make a marked contribution to the flow stress and hence hardness. Despite the importance of the ISE and its link to strain gradient plasticity (especially those variants based on GNDs) there have been very few independent measurements of GND density near indentations with which to validate the Nix and Gao theory. Indeed the understanding of flow fields around indentations is generally poor, with simple approximate models such as Johnson's cavity model often used.

Nye [5] gave an elegant geometrical link between lattice curvature and continuous distributions of GNDs described through the dislocation tensor. Nye's analysis was conducted under the assumption that there is no long range elastic strain field (i.e. stress field) present in the crystal. Removal of such an assumption introduces elastic strain gradients (i.e. stress gradients) into the formulation in addition to the lattice rotation gradients used by Nye, as has been discussed in detail by Kroner [START_REF] Kroner | Continuum theory of dislocations and self stresses[END_REF] .

Electron back scatter diffraction (EBSD) has been used by several groups to study localised plastic deformation. Early attempts used the contrast reduction and blurring of the EBSD patterns caused by defects within the diffracting volume as a measure of the extent of cold work [START_REF] Wilkinson | [END_REF]8], but with the advent of automated EBSD mapping interest has moved to correlating various measures of the misorientations found within the grains with imposed plastic deformation [9][10][START_REF] Lehockey | Electron Backscatter Diffraction in Materials Science[END_REF]. Adams and co-workers [START_REF] El-Dasher | [END_REF]13] were the first to attempt to use Nye's framework to recover GND densities from electron back scatter diffraction (EBSD) measurements of the spatial variation in crystal orientation. Their studies examined GND accumulation near grain boundaries in deformed Al bi-and poly-crystalline samples. Field et al [14] also used the method described by Adams and co-workers to determine the microscale GND distribution in high purity Al single crystals deformed in compression. Pantleon [15] has recently pointed out that earlier analyses [START_REF] El-Dasher | [END_REF]13] used too approximate a treatment for extracting the rotation terms and consequently lost some information that would further constrain the GND solution. Pantleon shows that five terms in Nye's dislocation tensor can be established directly from the EBSD measurements, and that the difference between two other terms can also be calculated and used as a further constraint on the GND solution. Pantleon applies his analysis to GND storage within a cold rolled Al polycrystal and demonstrates an improved lower bound solution is obtained when all six constraints are used. Wilkinson, Meaden and Dingley [16] have recently developed a new cross-correlation based analysis of EBSD patterns which allows rotations to be measured at much greater sensitivity of ~0.01° compared to 1° to 0.5° with 'conventional' Hough transform based analysis. This development removes much of the cumbersomeness in earlier work by Wilkinson [17,18] in which a large screen to sample distance was used to achieve the angular resolution. The increased sensitivity gives clear benefits in the recovery of GND density, which has been demonstrated by analysis of plastic deformation at the tip of a fatigue crack tip [19]. Further references concerning the EBSD technique can be found in the following reviews [20,21].

Here we report high resolution EBSD measurements of lattice curvature near nano-indents in Fe and determine from them GND distributions.

Experimental Method

The test sample was a Fe -0.01wt% C polycrystal obtained from Metal Crystals & Oxides Ltd. (Cambridgeshire, UK) with a grain size of ~ 50 µm. After polishing to 1 µm diamond paste the sample was electropolished in an electrolyte of 10% perchloric acid and 90% acetic acid at 19 volts and ~20°C for 45 seconds. Indents were made with a diamond Berkovich indenter tip using an MTS Nanoindenter XP. The indents were made under constant indentation strain rate of 0.05 s -1

, where the indentation strain rate is defined as indenter tip speed divided by indent depth. The orientation of the indented grain is illustrated by the wireframe inset in fig. 1a, and is described by the following Euler angles φ 1 =307°, Φ=32°, φ 2 =89°. EBSD measurements were made in a JEOL JSM 6500F FEG SEM operated at 20 keV, and a probe current of 2 nA. A peltier cooled CCD camera (Digiview 12) with 1300 x 1030 pixels and intensities digitised to 12 bit controlled by TSL OIM 1 software was used to capture the EBSD patterns. Measurements were made on the indented surface of the sample from which we obtained maps at a step size of 200 nm around the indents. Patterns were recorded at the full resolution of the camera for subsequent batchwise off-line analysis using CrossCourt 2.0 from BLG Productions Ltd 2 . Full details of the cross correlation analysis of EBSD patterns has been given in [16] including the effects of sub-region size, and noise on the patterns. The quality of the pattern shift measurements can be assessed using the height of the peak found in the cross correlation function. For a given sub-region the cross correlation peak heights are normalised to bring the peak height for the reference pattern (auto-correlation) to unity. If the pattern quality, that is the Kikuchi band contrast and sharpness, becomes reduced due to the presence of severe plastic deformation, or surface contamination the peak heights are reduced. Shadowing from surface topography also reduces the peak height for those sub-regions of the pattern falling within the shadowed part of the phosphor screen. The geometric mean of normalised peak heights for all sub-regions on a given pattern is calculated and used as a single parameter indicating how well the cross correlation pattern matching has performed. In the current study patterns giving a mean peak height value below a threshold value of 0.3 were rejected. This generally only happened in regions very close to the indent, where both shadowing and severe plastic deformation are present. Lattice rotations and elastic strains cause small shifts in the positions of zone axes and other features in the EBSD patterns obtained as the electron beam is scanned over the sample. These small pattern shifts are measured using automated image processing based on two dimensional cross-correlation analysis. The shifts at four or more sub-regions widely dispersed across the EBSD patterns are sufficient to directly calculate eight of the nine degrees of freedom contained in an arbitrary strain and rotation. The remaining degree of freedom is the hydrostatic dilatation of the lattice to which this method is insensitive, though it can be recovered if we, not unreasonably, assume that plane stress conditions are met within the diffraction region close to the sample surface. For the results presented here pattern shifts were measured at 20 sub-regions and least squares methods used to calculate the best fit displacement gradient tensor. The quality of this solution was assessed by calculating the angular difference between the shift expected at each of the sub-regions on the basis of the best fit solution and the actual measured shift. These differences were averaged over all 20 sub-regions to give the mean angular error, which provides a useful means of assessing how good the best fit solution is.

Analysis for GND Content

The analysis for GND content starts from the simple statement that we require continuity in the displacement field, so that integration of displacement (u) around any closed loop generates zero.

0 = = ∫ ∫ C C β dx du , (1)
where β β β β is the displacement gradient (or distortion) tensor, such that (2) For a given location in the crystal equations 1 and 2 must hold for any closed loop C, and for any surface patch S, from which it follows that the elastic and plastic parts of the distortion field must be related by ( ) ( )

j ij i dx β du = F o r P e e r
p e β curl β curl - = (3)
It is also significant to note that equation 1 is related to the Burger's circuit construction from which the net Burgers vector B of dislocation lines piercing the loop C, is obtained from

( ) ∫∫ ∫ - = - = S p C p β curl β ds dx B (4)
Remembering that Nye's dislocation tensor α α α α is defined by

∫∫ = S ds α B (5)
So the dislocation tensor can be found directly in terms of the elastic distortion field through ( )

e β curl = α (6)
Writing this explicitly in terms of derivatives of the elastic strain and rotation fields gives

                  ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ - ∂ ∂ +                   ∂ ∂ - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ =
ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε α (7)
High resolution EBSD measurements allow the rotation and strain field to be mapped over a two dimensional grid on the sample surface, however the variation with depth is not obtained. This means that terms containing 3 x ∂ ∂ are unknown leaving only three of the nine elements of the dislocation tensor fully determined. However, in many instances that we have looked at experimentally, we find that the rotation gradients are considerably larger than the elastic strain gradients, in which case we can either neglect the elastic strain gradients entirely or include those terms that are available experimentally and put the other terms to zero. In this case five of the nine elements of the dislocation tensor are determined directly, and additionally the difference between two elements (specifically α 11 -α 22 ) can also be found, as pointed out by Pantleon [15]. Pantleon's suggestion [15] that this difference between two elements in the Nye tensor can still be determined when elastic strain gradients are included appears to be in error, since the different symmetric and anti-symmetric nature of the strain and rotation tensors prevent the contributions from both of these sources being cancelled by taking the difference between terms in the Nye tensor. Recent work by Demir et al [22] use three dimensional EBSD measurements using serial sectioning in a FIB + SEM instrument to access the gradients in the third dimension. In principal this allows all nine elements in the Nye tensor to be determined under the assumption that elastic strain gradients can be neglected. The elements in Nye's dislocation tensor are related to the densities ρ s of different types of dislocations through

∑ = ρ = α max s s s j s i s ij t b 1 (8)
Where b s is the burgers vector, and t s the unit line direction of the s th dislocation type. Nye originally considered a set of nine dislocations with line directions and burgers vectors constrained to be along the cube axes of the crystal, in what Sun et al [13] refer to as the simple cubic (SC) deconstruction. In this case an unambiguous link between the nine elements of α and the densities of the different types of dislocations can be established. In our case this does not realistically reflect the crystallographic geometry of possible dislocation types in the BCC Fe crystal, however, such an approach may provide a lower bound estimate of the overall dislocation density. In our analysis for the SC deconstruction we will consider only rotation gradients and set terms that cannot be measured to zero. To better incorporate the crystallography of the Fe sample we can assume that the deformation is supported by movement and storage of dislocations with a/2<111> burgers vector slipping on {110} planes. We further assume that the dislocations are represented by either pure screw, or pure edge types. This results in 16 distinct dislocation types (ignoring positive and negative senses) including 4 screw types (with <111> line directions) and 12 edge types (with <112> line directions). As with the FCC case, which has been looked at by several other authors [13-15, 22, 23], there is considerably redundancy built into the system, so that there are many ways in which the measured values of α ij (and α 11 -α 22 ) can be supported by the available dislocations. This forces us to use optimisation routines which find dislocation densities that support the measured values of α ij and minimise some other related parameter. We write the problem with the dislocation densities represented by the column vector ρ ρ ρ ρ (16 by 1), and the measured rotation and elastic strain gradients making up elements of Nye's dislocation tensor also as a column vector Λ Λ Λ Λ (6 by 1 for analysis neglecting elastic strains, 3 by 1 when elastic strains are fully incorporated) as

Λ ρ A = (9)
The s th column of A contains the various terms in the tensor products of the Burger's vector and line direction of the s th dislocation type, as required by equations 8 and 7, and describe how an array of these dislocations would distort the crystal lattice. Two distinct possibilities arise in minimisations schemes that can be used to select a solution for the dislocations densities in equation (9). The mathematically simpler of these is the so called L 2 optimisation scheme which finds the solution giving the minimum sum of squares of the dislocation densities. This approach has been adopted by Field et al [14], and the solution for the equations can be written explicitly as

( ) Λ A A A ρ 1 - = T T (10)
Arsenlis and Parks [23] have explored these optimisation methods in detail for the FCC case in which 18 possible dislocation types are available. Assuming all nine terms in the Nye tensor are known they note that the null space of A gives nine linear combinations of dislocation densities ρ ρ ρ ρ s which together generate no contribution to the Nye tensor. Such combinations constitute statistically stored dislocations. In our case this situation is slightly different in that we have 16 dislocation types, but can measure only 6 of nine possible effects on the lattice distortion gradient. The null space of A is thus slightly larger, containing 10 combinations of the dislocation types that are either statistically stored or generate no distortion gradient that can be measured by the EBSD mapping on a single plane through the sample. The remaining six linear combinations of the dislocation types generate distortion gradients that are measured. The problem with the L 2 optimisation scheme is that there is no real physical argument leading us to expect a solution which minimises the sum of squares of the dislocation densities. We thus explore this scheme here mainly because others have used it and we wish to compare results with this analysis to those obtained by the L 1 optimisation scheme which we favour. In the L 1 optimisation we seek a solution to equation 9, that minimises the total dislocation line energy. All our dislocations have the same magnitude Burger's vector, so the only difference in their energies is generated by the ratio of energies for edge and screw dislocations Other groups have sought to minimise the total dislocation density [START_REF] El-Dasher | [END_REF]13,22] or equivalently line length. This is very similar to the approach taken here but in effect the ratio of energies of edge and screw dislocations are taken as unity rather than ~1.5 which is more typical of a metal (ν=1/3). Various numerical schemes are available for solving this type of problem including the classical simplex algorithm. Here we use the optimisation algorithms implemented in MatLab

3

. Again the solution found represents one of many possible solutions that satisfy the minimisation, and equation 9.

Results

The crystal orientation and the geometry of the indent with respect to the reference axes used are given in figure 1, along with a typical EBSD pattern. Figure 2 compares parameters used to assess the quality of data obtained using the conventional Hough transform based analysis of the EBSD patterns with their closest equivalents available from the newer cross-correlation based analysis. The image quality map from OIM and the geometric mean peak height map from CrossCourt indicate where pattern contrast is reduced, and where shadowing arises from the topography of the indenter impression on the scanned surface. The variation at the indent edge is much more abrupt for the OIM image quality map (fig 2a). By comparison the peak height map (fig 2c) shows a smoother transition at the indenter edge but continues to give some variation across the flat surface, with lower values generally recorded within the deformed regions close to the indent. The mean error in interplanar angles from the Hough based analysis (fig 2b) also shows a step change at the edge of the indent with an average value of ~0.5°, (or ~8x10 -3 rads) obtained outside the indent. The mean angular error in the cross-correlation based shift measurements are ~100 times smaller; the average value outside the indent being ~0.005°, or 9x10 -5 rads. In the region marked A in fig. 1a the peak height is somewhat reduced while the mean angular error is slightly raised compared to other regions around the indent. This small effect is caused by surface topology due a slight pile-up of material in this region. It should be noted that previous analysis of the cross-correlation method indicate that the errors are ~5 times larger on some components of the strain and rotation tensors than on those with the smallest measurement errors [19,24]. Thus although the mean angular error gives an indication of how good the best fit solution is the measurement noise will be greater on some components of strain and rotation than others.

Rotation Fields

One expects the lower error indicated in the cross correlation based analysis to be linked to less noise and greater detail being visible in the maps of the rotation field. This is indeed the case as is clear in figure 3 which shows the rotation field obtained for the Hough transform based analysis (fig 3 a to c) and the cross correlation based analyses (fig 3 d to f) for the same 500 nm deep indent shown in figures 1 and 2. Data are only shown where the mean error in interplanar angles is 1° or less in the Hough transform based analysis, and where geometric mean of the cross correlation peak heights was found to be 0.3 or greater in the cross correlation based analysis, so that regions near and within the indent where the patterns were of poor quality are excluded. The rotation has been decomposed into 3 components about the orthogonal reference axes (see fig 1) using the assumption of small angles and linear superposition (note most EBSD measurements involve larger angles which requires the more complicated multiplicative decomposition). Near face A we expect the plastic flow away from the indenter face to cause lattice rotations that are predominantly about the x 1 axis (ie the w 23 The rotations components generally take values considerably larger than the mean angular error in pattern shifts obtained with the cross correlation analysis. In contrast, the errors associated with the Hough transform based analysis are of comparable magnitude to the rotation fields being measured, and so although similar variations in rotation are seen compared to the cross correlation based analysis the noise level is significantly higher. Noise on the orientation measurements affects the uncertainty in the misorientation axis and becomes very marked as the misorientation angle becomes comparable and indeed smaller than the orientation sensitivity [25,26] . Bate et al [26] derived the following equation for the uncertainty in the misorientation axis φ, in terms of the misorientation angle θ and the uncertainty in the misorientation angle δ φ = arctan(δ/θ) [START_REF] El-Dasher | [END_REF].

Maps of this uncertainty in the misorientation axis is shown in figure 4. For the Hough transform generated data (fig 4a) the mean error in the interplanar angles was used as an estimate of δ, and given the small size of the misorientations present in the data very large uncertainty is generated in the misorientation axis information. It is only in the most highly misoriented regions that there is any real confidence in the misorientation axis. In contrast the cross correlation generated data (fig 4b) retain higher confidence in the misorientation axis across the whole field because of the much smaller error in the misorientation angle δ (taken to be the mean angular error in the pattern shifts).

Elastic Strain Fields

The cross correlation analysis has sufficient sensitivity for the variation in elastic strains to be mapped in addition to the lattice rotation. Maps of selected components of the elastic strain tensor for this indent are given in figure 5, again data is only shown if the geometric mean of cross correlation peak heights is 0.3 or above. The strains are generally smaller than for the rotations and so these maps appear a little noisier than those in figure 3, however, systematic elastic strain variations are clearly resolved by the method. The elastic strains are referred to polar coordinates with the origin marked on the figures at the centre of the indent. The radial strains are seen to be generally compressive, while the hoop strains are generally tensile which is in line with expectations. The triangular pyramid shape of the indenter imparts some angular variation of these strains not captured by axisymmetric models such as Johnson cavity model. These angular variations of strains is further promoted by the crystallographic anisotropy of plastic slip relieving stresses in some directions more effectively than others. The in plane ε rθ shear strains not accounted for in axisymmetric models are also present at significant levels, though are limited to small regions close to the indent and fall of quickly to small values at larger radii. The other elastic strain components, that is normal, and shear strains containing out of plane displacements are significantly smaller and are not shown here.

GND Content

Various routes for constructing maps of the GNDs density distributions near the indent using the Nye/Kroner analysis have been explored. Summary results showing the total GND density fields are compared in figure 6. The maps in figures 6a and 6b were constructed from Hough based analysis of the lattice rotation field. It is difficult to discern much obvious structure to the storage of GND in the plastic zone near the indent. A simple lower bound estimate of the expected noise on these GND measurements, can be made from the uncertainty in the orientation measurements δ = 8x10 -3 rads and the step size used in the EBSD map λ of 200 nm using the simple expression

λ δ ρ b = ∆ (13)
Which gives a large uncertainty in the GND density of ~2x10 14 m -2

, which masks the underlying distribution. GND distributions generated from the cross correlation based analysis of the EBSD patterns (figure 6 c to h) exhibit a reduced noise level, estimated to be ~2x10 12 m -2 and allow the plastic zone to be revealed. This simple estimate ignores the influence of misorientation axis and the manner in which errors are compounded by the addition of densities on many slip systems. Analysis of EBSD maps on well annealed samples away from any indent indicate a noise level of ~10 13 m -2 for the total GND content under the BCC deconstruction. Comparing figures 6 (c) and (d) shows that the simple SC deconstruction generates a significantly lower estimate of the total GND density than the BCC deconstruction. Curvatures are less efficiently accommodated in the BCC deconstruction because the true crystallography leads to a system of dislocations with Burgers vectors (and indeed line directions) that are not mutually orthogonal. The two optimisation methods are compared in figures 6 (d) and (e), where it is seen that in this instance the L 2 scheme which has less physical basis leads to a very similar pattern for the GND distribution and generates densities that are marginally higher than for the L 1 scheme which we prefer. As expected from the smaller magnitudes of the elastic strains, inclusion of elastic strain gradients as corrections to the curvatures calculated from the rotation field only leads to small changes in the estimate of the GND density. This can be seen by comparing figures 6 (d) and (f) though here not all terms required to fully correct for the strain gradients are available. If we limit the analysis to just three terms in the Nye tensor then it is again seen that inclusion of the elastic strain gradient has only a small effect (cf figures 6 g and h). A much more significant effect on the GND estimate is that using only three terms in the Nye tensor leads to significantly lower values for the GND density (cf figures 6 f and g, and figures 6 d and h). Our preferred method is to use the six available derivatives of the lattice rotations corrected for the available elastic strain gradients within an L 1 optimisation to minimise the total line energy (ie fig 6f). We also note that having some of the elastic strain fields measured allows us to assess whether or not the neglect of some elastic strain gradients is likely to lead to large errors in determining the GND content. Figure 7 shows the distribution of each of the 16 individual GND types making up the BCC deconstruction. It should be noted that there is a non-uniqueness issue within this analysis so that this is one possible solution that satisfies the measured lattice curvatures and minimises the total dislocation line energy. The fact that regions of higher density for particular GND types are spatially correlated gives some support to the idea that this particular solution is a reasonable one. Furthermore the dislocation densities appear to be of reasonable magnitudes. In figure 7 negative dislocation densities simply mean that the sense of the dislocation is the opposite of that indicated by the 'clockface' found toward the top left of each map. This 'clockface' indicates the projection of the Burger's vector onto the sample surface plane. A Burgers vector fully in the surface plane reaches to the edge of the circle while in other cases the component in the surface plane determines the length, while colour (on-line version) is used to indicate whether the out of plane component is into (red) or out-of (green) the sample surface. A similar approach is taken for the line direction shown black if out-of the sample and blue if into the sample. It is noticeable that the lower line energy of the screw dislocations generally leads to them being preset in greater densities (fig 7 a to d) than the edge types (fig 7 e to p). In many cases the regions of high density for a particular dislocation type are in bands that radiate out from the indent. In some instances (eg fig 7g) there are several parallel bands extending out from different parts of the indent. Figure 7 suggests that some dislocation types are considerably less populated than others, notably the edge dislocation in fig 7k andfig 7m, and of the screw dislocations those in fig 7d. There are several possibly reasons for the lower populations of some types of GND in the optimal solution: (i) these dislocations are simply not present, (ii) the geometry of the dislocations leads to curvatures that are not readily detected, or, (iii) these dislocations may be present but are not reflected by the optimisation to minimise line energy. Examining the magnitudes of terms in the columns of matrix A helps to rule out possibility (ii) as we find that in this particular crystal orientation all the dislocation types should generate measureable curvatures of similar size. Indeed the norms of columns making up A have an average of 0.19 and a standard deviation of 0.04.

Discussion

The cross correlation analysis of the EBSD patterns clearly achieves higher sensitivity than the more standard Hough transform based analysis. Visual comparison of the lattice rotation fields obtained from the Hough-based and cross correlation-based analysis (fig 3) tend to suggest that the Hough-based analysis is able to pick up the underlying trends despite the significantly greater noise level. However, the uncertainty in the rotation axis (fig 4) is marked and produces a large uncertainty in determining any information concerning the Burgers vector and line direction of the dislocations giving rise to the rotation. The need for lattice rotation (and elastic strain) gradients in the analysis for GND content exacerbates the noise problem since (numerical) differentiation tends to always leads to increased noise. With the step size used here a simple analysis leads to a lower bound estimate of the noise on the GND density of ~10 14 m -2 for the Hough-based analysis. This is already well into the densities expected for quite extensively cold worked metals so that only the very extreme GND densities can be reliably measured. At first glance equation 13 suggests that increasing the step size for the EBSD map would lead to increased sensitivity in the GND density measurement, however, increasing the step size also increases the likelihood of there being combinations of dislocations between the measurement points that make no (or reduced) contribution to the lattice curvature. In other words keeping the step size small has the attraction that a larger fraction of the dislocation content can be thought of as GNDs which we can detect, rather than statistically stored dislocations which are missed in the analysis. Attempts to reduce noise by using more than two measurement points to define the orientation gradient would , for a given step size, also increase the effective length scale over which the GND content was assessed and may omit real features present at shorter lengths. Such step size dependence of more empirical measures of deformation from EBSD determined misorientation maps have been explored by LeHockey [START_REF] Lehockey | Electron Backscatter Diffraction in Materials Science[END_REF] and Kamaya et al [9]. In comparison at ~10 12 m -2 the lower bound estimate of the noise level on GND content derived from the cross correlation analysis of EBSD patterns is such that dislocation structures in deformed or even annealed metals can be assessed, as indeed can strains and defect content in semiconductor materials and structures. The range of applications of this cross correlation based analysis of EBSD patterns is large and so it is important to generate parameters that allow the quality of patterns and data extracted from them to be assessed in an automated and quantitative way so that the rapidly growing size and number of linescans and maps that are generated can be used with confidence. Two such parameters, namely the geometric mean of the normalised cross correlation peak heights and the mean angular error in the measured pattern shifts usefully perform this function. Figure 2 demonstrates that they can both be used to distinguish between regions in a map where the data is sound and where it is of poorer quality, in this case inside the indent. Systematically mapping dislocation densities in the range 10 12 to 10 16 m -2 over areas of the size studied here or even larger, would be too labour intensive to undertake using TEM and would require painstaking or possibly prohibitive sample preparation. Conventional x-ray and neutron diffraction methods can probe dislocation densities in this range using peak broadening analysis; however, this provides only volume averaged information. The EBSD method provides a means of probing dislocation arrangements at the mesoscopic and microscopic length scales and so is an important bridge between analysis of individual defects in the TEM and volume averaged densities available through x-ray and neutron probes. It is of course at these mesoscopic length scales that many important phenomena occur during plastic deformation as a result of complex interactions between dislocations and microstructural features. Some synchrotron based x-ray diffraction techniques [27][28][29] do provide information that is more directly comparable to that obtained with EBSD. These approaches have the advantage that the full three dimensional nature of the problem can be addressed, but have the disadvantages of lower spatial resolution and the inconvenience and cost associated with large centralised facilities with limited access. Vaudin et al [30] have recently used EBSD patterns and CrossCourt software to measure deformation approaching a wedge indent in Si. To ease direct comparison with micro-Raman measurements from the same indents they present the EBSD results for the variation of in plane biaxial stress, which is shown to follow a 1/r 2 dependence, with r being the distance from the central line of the indent. This is the variation expected from a simple comparison with the elastic strains and stresses close to an expanding cylindrical cavity. The EBSD results also agree extremely well with micro-Raman measurements when suitable exciting illumination is selected to minimise the penetration into the sample depth. In our case the analogous elastic analysis is of a spherical cavity a leads to a 1/r 3 variation of elastic strains. However, there is markedly more plasticity near the indent in Fe compared to Si and so the resulting elastics strains are smaller and do not show as simple a trend. Kysar et al [31] have used Hough-based EBSD to examine sub-surface deformation in cross sections through relatively deep (several hundred microns) wedge indents in Cu and Al single crystals. Lattice rotations below these indents were significantly larger, at up to 20°, than those seen here at the free surface for much smaller indents. Plane strain conditions were assumed in a GND analysis of the lattice curvatures leading to a reduced set of dislocation types the densities of which are found using the less physical L 2 minimisation scheme. Although the symmetry of the crystals and indenter do imply plane strain conditions at a macro-scale at the mesoscopic scale probed by EBSD this may not be the case. It may be of interest to use a more general analysis such as that presented here and use this to examine the extent to which constraints on the deformation expected from a macroscopic viewpoint are in fact followed at the more local level within the dislocation sub-structure. Hough-based EBSD studies of sub-surface deformation for indents in Cu have also been made by Demir et al [22]. They studied conical indents at depths from 450 nm to 1230 nm, using repeated sectioning in a FIB-SEM instrument to probe the 3 dimensional structure. Again large rotations are seen below the indents making the Hough-based analysis possible, though they note a noise level of ~10 densities in accord with our estimation here. Rester et al [32] have also used Hough-based EBSD to map orientation changes below cube-corner indents in Cu and note that the self-similar pattern in the misorientations breaks down for very shallow indents below ~300 Nm deep, and indeed for very large indents to depths greater than ~30 µm. The Nix-Gao theory [4] for the ISE assumes that the GNDs are stored in a hemispherical region below the indenter contact in a uniform distribution. Furthermore for convenience the GNDs are considered to be arranged as prismatic loops with Burgers vectors perpendicular to the surface which are rather efficient at accommodating permanent impression made by the indenter. For a Berkovich indent of 500 nm depth in Fe the Nix-Gao model gives a GND density of ~1.5x10 15 m -2 and densities close to this were determined from the EBSD measurements, and may be even higher below the indented surface. In the experiment it is clear that the GNDs are far from uniformly distributed despite the availability of many slip systems in BCC structure. In some place the high GND densities extend a considerable distance beyond the indent impression, and the predominant burgers vector varies around the indent in a complex way. The detailed interaction of the triangular indenter and motion and accumulation of dislocations in the crystal is not captured by the Nix-Gao model which necessarily simplifies the geometry to generate a tractable problem. Crystal plasticity finite element analysis by Liang et al [START_REF] Liang | International Conference on Heterogeneous Material Mechanics[END_REF] has been used to simulate the deformation around the indent studied here. The ISE is captured in the simulation which includes hardening from GNDs which are allowed to evolve (increase and decrease) in response to changes in the plastic strain gradient. Lattice rotations are predicted by the model and are in qualitative agreement with experiment. The magnitudes of the rotations are slightly over predicted but this may be due to omission of the unloading phase from the simulations.

Conclusions

We have used a cross correlation based analysis of EBSD patterns to map the variation of the elastic strain tensor (and hence stress tensor) and small lattice rotations in the vicinity of a 500 nm deep Berkovich indent in Fe. Two measures of the data quality have been developed which successfully identify regions where shadowing and severely blurred patterns prevent sensible analysis from being undertaken. These are the geometric mean of the normalised cross correlation peak heights and the mean angular error in the measured pattern shifts. Lattice rotations are significantly larger than elastic strains. The radial elastic strains are generally compressive and the hoop strains tensile as expected, though the variation is strongly influenced by the anisotropy in the plastic flow. Direct comparison of maps from the cross correlation based analysis with those from the conventional Hough-based analysis clearly demonstrate the improved sensitivity of the cross correlation technique. Analyses for the GND content was undertaken assuming dislocation types representative of both simple cubic and body centred cubic crystals and using as inputs various components of the gradients of the lattice rotations and elastic strains. Our favoured approach is to seek a combination of GND types which support the 6 measureable (of a possible 9) gradients of the lattice rotations after correction for the 10 measureable elastic strain gradients, and minimise the total GND line energy using an L 1 optimisation method. The corrections for elastic strain gradients neglect 3 terms that cannot be measured as gradients along the surface normal are not accessible. The solution thus forms a lower bound estimate of the GND density. The ±10 -4 sensitivity of the cross correlation analysis corresponds to a lower bound estimate of ±10 12 m -2 for the noise on the GND density for a 200 nm step size. GND densities as high as 10 15 m -2 were found near the indent. The lower line energy of the screw dislocations makes them more favoured in the solution, and regions were specific dislocation types are located are observed in the solution.

  component, (fig 3 c and f), and that is seen to be the case. The other two components, w 31 about the x 2 axis (fig 3 b and e ) , and w 12 about the sample surface normal x 3 (fig 3 a and d), are both generally smaller than w 23 and vary considerably with position along face A, even changing sign from one end to the other. Compared to near face A, the w 23 rotations are smaller and of opposite sign near faces B and C.Along the plane bisecting face A parallel to the x 2 axis the indenter has mirror symmetry, however, this symmetry is not maintained in the distribution of w 23 rotations. This loss of the indenter symmetry is even more marked for the w 12 and w 31 components of rotation, and is presumably due to the specific alignment of the available Burgers vectors and slip planes in the indented crystal. The four possible a<111>/2 burgers vectors are shown in fig1, and it is clear that these do not have mirror symmetry through face A.

Figure 1 :

 1 Figure 1: (a) secondary electron image of 500 nm deep Berkovich indent in Fe with isets showing the crystal orientation and the trace of <111> slip directions in the surface plane. (b) An EBSD pattern from the crystal with sub-regions used for cross correlation analysis marked.

Figure 2 :

 2 Figure 2: maps showing EBSD data quality. (a) image quality index based on height of peaks found in Hough transform, (b) mean error in interplanar angles from Hough based analysis (in degrees), (c) geometric mean of normalised cross correlation peak heights, and (d) mean angular error in pattern shifts determined in cross correlation based analysis (in degrees).

Figure 3 :

 3 Figure 3: lattice rotation fields measured by Hough based analysis (a to c) and cross correlation based analysis (d to f). Rotation angles are in radians.

Figure 4 :

 4 Figure 4: Estimated error in alignment of misorientation axis (in degrees) for (a) Hough based analysis and (b) cross correlation based analysis.

Figure 5 :

 5 Figure 5: Selected elastic strain components given in cylindrical polar co-ordinates referred to an origin marked at centre of the indent. Other components are also determined but are markedly smaller.

Figure 6 :

 6 Figure 6: Estimates of GND density distribution near indent. [Colourscale shows log(ρ) with ρ in lines/m -2 .] (a) Hough based rotation fields, strains not measured, SC deconstruction of 6 terms in α. (b) Hough based rotation fields, strains not measured, BCC deconstruction of 6 terms in α. (c) cross correlation based rotation fields, strains neglected, SC deconstruction of 6 terms in α. (d) cross correlation based rotation fields, strains neglected, BCC deconstruction of 6 terms in α, L 1 optimisation. (e) cross correlation based rotation fields, strains neglected, BCC deconstruction of 6 terms in α, L 2 optimisation. (f) cross correlation based rotation fields, corrected for measured elastic strain gradients, BCC deconstruction of 6 terms in α, L 1 optimisation. (g) cross correlation based rotation and strain fields, BCC deconstruction of 3 terms in α, L 2 optimisation. (h) cross correlation based rotation fields, BCC deconstruction of 3 terms in α, L 2 optimisation.

Figure 7 :

 7 Figure 7: One possible solution for densities of the 16 different GND types assumed in the BCC deconstruction of cross correlation based rotation fields, corrected for measured elastic strain gradients to give 6 terms in α. L 1 optimisation to give a minimal line energy solution was used. The 'clockfaces' show the projection of the Burger's vector and line directions onto the sample surface plane. Burgers vectors are red or green to indicate a component along x 3 that is into or out-of the surface plane. Line directions are shown black if out-of the sample and blue if into the sample. (a) to (d) show the four screw dislocations while (e) to (p) the 12 edge types.

  Figure 1

Figure 5

 5 Figure 2

  Figure 6

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

MatLab, The Mathworks Inc (www.mathworks.com)

Acknowledgements

Thanks to (i) Dr Mark Vaudin (NIST, Gaithersburg) for discussion of EBSD results on indents in Si prior to publication and initial MatLab script for plotting maps from CrossCourt output, (ii) Dr Graham Meaden (BLG Productions Ltd) and Prof David Dingley (University of Bristol) for continuing discussion of all aspects of EBSD (iii) Prof Fionn Dunne (University of Oxford) and Dr Hong Liang (University of Oxford) for insightful discussions on their crystal plasticity simulations of indentation, and (iv) Ben Britton (University of Oxford) for discussions on GND analysis and extensive input into Matlab coding.