## Abstract

Recent experiments indicate that the first yield point of micron-sized metals exhibits significant statistical scatter as well as strong dependence on the specimen size. In this work, molecular dynamics (MD) simulations are carried out to investigate the onset of shear deformation in a small block of material, using an embedded atom potential for the intermetallic Ni_{3}Al alloy. Incipient plasticity in the form of homogeneous dislocation generation is observed to occur at atomic sites with interatomic displacements approaching one-half of the Shockley partial Burgers vector. From the distribution function of the interatomic displacements observed in the MD simulations, the probability of a general material volume surviving under given loading conditions is predicted. The survival probability is then calculated for several situations, including homogeneous deformation and nanoindentation, to predict the critical load for incipient plasticity to occur in these situations. The predicted results are compared to micro-pillar compression and nanoindentation experiments on Ni_{3}Al available in the literature.

## 1. Introduction

Recent rapid advancements in micro-machinery technologies call for an urgent need to understand the mechanical behaviour of materials of dimensions in the sub-micron regime. It is by now well known that a wide range of crystalline materials exhibit the so-called size effects of strength. These refer to an increase of the flow stress at some finite strain, in terms of hardness obtained from nanoindentation measurements (Pethica *et al*. 1983; Doerner & Nix 1986; Stelmashenko *et al*. 1993; Ma & Clarke 1995; McElhaney *et al*. 1998; Lim & Chaudhri 1999; Liu & Ngan 2001), or torsional strength obtained from twisting tests on small whiskers (Fleck *et al*. 1994). An explanation, which has gained rather widespread acceptance, is that in these deformation configurations, strain gradients exist (Fleck *et al*. 1994), and these become more significant as the deformation size decreases, so that the working hardening effects imposed by the geometrically necessary dislocations accommodating the strain gradients are stronger at smaller deformation sizes (Nix & Gao 1998).

On the other hand, the onset of permanent deformation from a hitherto virgin state, i.e. the so-called incipient plasticity process, has also been found to exhibit significant size effects. Early tensile experiments on annealed metallic whiskers revealed that these small specimens exhibited sharp yield points which would not be present in bulk samples (e.g. Smallman & Bishop 1995). Also, the yield points were observed to occur at very high applied stresses. Similar observations were also made recently from compression tests on micro-pillars fabricated by the focused-ion-beam technique (Uchic *et al*. 2004). These experiments showed that micro-pillars of Ni and Ni_{3}Al exhibited a marked increase in the initial yield strength as the specimen size decreased. In uniform tension or compression tests, significant strain gradients would not be present and so geometrically necessary dislocations should not play any significant role in these situations. The scarcity of source dislocations is thought to be the reason for the much higher yield points in these small samples (Smallman & Bishop 1995), but the gradual reduction in the yield point as size increases as observed by Uchic *et al*. (2004) still lacks a satisfactory explanation.

Another aspect of the incipient plasticity of micron-sized crystals is that its occurrence usually exhibits significant statistical scattering. Incipient plasticity is easily observable in nanoindentation experiments on annealed crystalline materials as a load drop or displacement jump, and such experiments on annealed Fe–3%Si (Gerberich *et al*. 1995), Ni_{3}Al (Chiu & Ngan 2002; Wo & Ngan 2004; Wo *et al*. 2005), and Al (Feng 2001) revealed that, by holding the applied load initially in the elastic range, yielding may occur subsequently and suddenly after a certain waiting time, implying that these materials cannot indefinitely sustain GPa-level stresses in the elastic range. The waiting time for incipient plasticity to occur, however, was not constant even at a fixed temperature or load, but exhibited a statistical distribution (Gerberich *et al*. 1995; Feng 2001; Chiu & Ngan 2002; Wo *et al*. 2005). When the holding load was high, the distribution function of the waiting time was found to follow exponential decay, so that the chance of longer waiting time, or sample lifespan if plasticity is regarded as failure, would be smaller. When the holding load was low, the probability function of the waiting time was found to be a peaked function with a modal value decreasing at lower holding loads (Wo *et al*. 2005). Even in constant-load-ramp nanoindentation experiments, the transition from elasticity to plasticity does not occur at a fixed load at a given temperature, but instead exhibits significant statistical scatter (Schuh & Lund 2004; Wo & Ngan 2004; Wo *et al*. 2005).

In terms of both the size dependence and the stochastic nature, the yield strength of small samples is rather similar to the fracture strength of bulk ceramics, which is well known to exhibit significant statistical scatter as well as sharp dependence on the sample size. The strength of bulk ceramics is governed by the spatial distribution of pre-existing cracks, and, more importantly, the chance of encountering them in a given material volume and stress state. A well-established methodology to design against the size dependence and the stochastic nature in bulk ceramics is the Weibull statistics (Weibull 1951). It would thus appear that a model similar to the Weibull statistics can be developed to describe the size dependence and the stochastic nature of the incipient plasticity in micron-sized crystals.

The Weibull statistical model developed for bulk ceramics, however, is purely phenomenological. In the present work, rather than using a similar phenomenological approach, we use a scaling approach to predict the strengths of larger, micro-sized volumes, from molecular dynamics (MD) data on smaller, nanometre-sized volumes. This is based on the assumption that the probabilistic distributions of atomic properties, such as relative displacements, energies, and so on, are the same for small and large samples, as long as plasticity has not occurred so that the atomic arrangement is still perfectly crystalline.

## 2. MD simulations of dislocation generation

MD simulations were carried out: (i) to study how homogeneous dislocation nucleation occurs in Ni_{3}Al under intense shear and (ii) to investigate whether scaling relations can be obtained to predict nucleation in much larger material volumes. The Voter–Chen embedded atom method (EAM) potential for Ni_{3}Al (Voter & Chen 1987) was used in the MD simulations. The simulation block had 36×36×24 atoms, measuring 4.5×7.4×10.5 nm, along the three crystallographic directions of , [111] and , respectively. Periodic boundary conditions were used along the and directions, and the (111) sides of the block were simulated as surfaces bounding the block. A simple shear strain was applied by keeping the bottom (111) layer stationary and moving the topmost (111) layer at a constant velocity along the direction throughout the simulation (Horstemeyer *et al*. 2001), and setting an initial velocity gradient between these two layers. The shear strain rate of the simulation was set to be 5×10^{8} s^{−1}. The MD code used was based on the leap-frog formulation of the Verlet algorithm (Allen & Tildesley 1987), with the initial velocities set up according to the Maxwell distribution. The time step for the simulation was 6 fs, so that a vibrational period was split into about 10 steps to ensure that velocities and accelerations were constant within a time step. At this time step, the typical fluctuations in the total energy during simulation were smaller than 1 part per 5000 over 20 time steps. The temperature was controlled at different values ranging from 20 to 600 K during the simulations using a Woodcock thermostat (Woodcock 1971). Incipient plasticity was detected by looking for sudden drops in the potential energy and stress, the latter was calculated using(2.1)where *τ*_{αβ} is an off-diagonal element of the stress tensor, *V* is the volume of the simulation cell, *m*_{i} and *p*_{iα} are the mass and the *α*-component of the momentum of the *i*th atom, and *r*_{ijα} and *f*_{ijβ} are the *α*-component of the distance and the *β*-component of the force between atoms *i* and *j*, respectively (Allen & Tildesley 1987). Figure 1 shows the potential energy as well as the shear stress component along the applied shear strain direction in a typical constant-rate shearing simulation, from which clear drops in the potential energy and stress can be observed at simulation steps close to 33 560, corresponding to an instantaneous strain rate of about 0.1. The sudden drops in potential energy and the shear stress correspond to a major relaxation event accompanying incipient plasticity in the crystal.

To scrutinize the actual incipient process, atoms having high magnitudes of maximum relative displacement with respect to their neighbours are traced at simulation steps just before and after the sudden relaxation event. The maximum relative displacement *x*_{i} of each atom *i* is calculated as(2.2)where *j* includes all the nearest neighbours of *i*, is the relative position vector of *i* and *j* and is the relative position vector at absolute zero temperature (Ngan *et al*. 2004). The more energetic atoms in the simulation cell can be revealed by displaying only those atoms with *x*_{i} higher than a certain cut-off value. Figure 2*a*–*d* shows the positions of the atoms with *x*_{i} larger than 17% of the lattice constant *a*_{o} at different snapshots during the shearing process. The cut-off value of 0.17*a*_{o} chosen here is quite close to one-half of the magnitude of the 1/6〈112〉 Shockley partial vector, which is 0.2*a*_{o}, and so the atoms displayed in figure 2 are among the most severely displaced in the simulation block. At an applied strain significantly below the critical value, energetic atomic clusters with large *x*_{i} appear randomly in the simulation block, as shown in figure 2*a*. These clusters, which are referred to as ‘hot spots’ hereafter, will not last for more than 10 simulation time steps, corresponding to 60 fs, which is the order of magnitude of the atomic vibration cycle. At a low strain, new hot spots appear at new sites after old ones are annihilated, and the hot spots have rather diffused sizes depending on the critical relative displacement used to define or display them. These hot spots are, therefore, likely to be momentary extrema in the lattice vibration and they are certainly not atomic defects. However, at a high enough strain close to the critical value, one or two hot spots may evolve during their limited lifetime of existence into a small dislocation loop as shown in figure 2*b*,*c*, and this quickly expanded to cover the whole slip plane in the simulation block as in figure 2*d*. The dislocation loop was identified to be a faulted loop with Burgers vector . In the L1_{2} crystal structure, the complex stacking fault bounded by the Shockley partial is a high-energy fault. The eventual wide expansion of the loop shown in figure 2*d* is therefore not very realistic, as this is likely to be due to the attractive forces imposed by the periodic boundary conditions. However, the early stages of the nucleation from figure 2*a*–*c*, comprising the evolution from a hot spot into a small loop, should be realistic enough, since cell boundary effects should be small on the hot spot and the small nucleated loop. In reality, once a small Shockley partial loop is nucleated, further faulting may occur which may be followed by cross slip to form a Frank–Read source for a subsequent, observable strain burst. In the following, we assume that the initial faulting process as shown in figure 2*a*–*c* is the critical process for a strain burst. This assumption is based on the recent experimental observation by the present authors (Wo *et al*. 2005) that at high enough nanoindentation loads, the statistical variation of the incipient plasticity process in Ni_{3}Al obeys homogeneous nucleation kinetics, but the measured activation volume is atomic rather than the volume expected from a dislocation loop. The atomic activation volume is likely to correspond to the hot spots in the present simulations, and hence the evolution of a hot spot into a small loop, instead of the subsequent processes, should be the activation process for the entire strain burst phenomenon observed in Ni_{3}Al.

Having established that the hot spots are precursors to homogeneous dislocation generation, we next analyse the probability of their formation. One very important discovery from the present simulations is that at any simulation step before the onset of defect generation, the fraction of atoms having a maximum relative displacement *x*_{i} higher than a given value *x* is well describable by the Weibull form,(2.3)where *x*_{o} and *m* are constants at a given simulation step. In our simulations in which the strain rate was constant, the strain rate used was sufficiently slow before the generation of dislocations, and so the shearing process can be regarded as quasi-static. Under this assumption, the simulation step is directly proportional to strain. The value of the normalization constant *x*_{o} is found to change significantly at different simulation steps or applied strains, as well as at different simulated temperatures, but the modulus *m* is found to remain within a narrow range of about 3.7±0.3 in all the simulations from 20 to 600 K. Figure 3 shows some examples of the excellent fit of equation (2.3) to a typical set of simulated results at different simulation steps and temperatures under constant-rate shearing. To see if the simulation block size has any significant effect on the validity of equation (2.3), the simulation at 300 K was repeated using a larger block of 46×36 ×24 atoms, which was approximately 30% larger than that used in the main simulations. It was found that the displacement profile of the atoms can also be very accurately described by equation (2.3), and the changes in the values of *m* and *x*_{o} were only within 14%. It is thus apparent that the block size used in the main simulations was already large enough to achieve reasonable convergence.

From equation (2.3), the mean relative displacement of the atoms, 〈*x*〉, is given by(2.4)where . Hence *x*_{o} is directly proportional to the mean relative displacement 〈*x*〉 if *m* is approximately constant. For example, if *m*=3.7, 〈*x*〉=0.9025 *x*_{o}. 〈*x*〉 can be calculated directly from the simulated configuration at each time step as the mean value of *x*_{i} in equation (2.2). Figure 4 shows the variation of 〈*x*〉 at different simulation steps and temperatures, from which 〈*x*〉 can be concluded to vary linearly with the applied strain, while the slope and intercept of this variation are functions of temperature. With equation (2.4), the results in figure 4 suggest the following phenomenological relation for *x*_{o},(2.5)where *γ* is the simulated shear strain along the system, and *A*(*T*) and *B*(*T*) are functions of temperature *T* only. Curves fitting the results in figure 4 yield the following approximate expressions for *A* and *B*:(2.5a)(2.5b)Physically, equation (2.5) suggests that the mean relative displacement between atoms has two components: (i) the *Aγ* component representing the biased relative displacements along the straining direction due to the applied strain and (ii) the *B* component representing the more isotropic, random relative displacements due to thermal agitation. The parameter *A* can therefore be expected to be proportional to the shear modulus of the material along the straining direction, and its slight but decreasing trend with temperature is likely to be associated with modulus softening at high temperatures. The increasing trend of parameter *B* with *T* can be interpreted as a reflection of increasing random agitation at higher temperatures.

## 3. Scaling model for elastic instability

Although the strain rate of 5×10^{8} s^{−1} used to obtain equation (2.3) is unrealistically high compared to real situations, it is nevertheless slow enough to allow quasi-equilibrium to be reached between successive rises in strain during the MD simulation. In other words, the probability function in equation (2.3) can be expected to be applicable to other strain rates, including perhaps real-time rates. In real-time experiments, the shear strain *γ* encountered would be much lower than that used in the MD simulations, and from equation (2.3), this means a much lower probability *P* for a given displacement *x*. Thus, the use of equation (2.3) to predict real-time situations is essentially an extrapolation scheme which bridges a big gap between the MD and the real-time regimes. Apart from the various assumptions behind the MD simulations, the validity of this extrapolation scheme remains a basic assumption in the present analysis.

As similar to the Weibull statistics for bulk ceramics, we focus on a probabilistic rather than a deterministic description for incipient plasticity. Let *F* be the probability of a sample surviving, i.e. remaining elastic, under a given applied stress state. Dividing the material volume into elemental volumes d*V*_{i}, the probability that the entire sample survives is , where *F*_{i} is the probability that d*V*_{i} at site *i* can survive. Hence,(3.1)Imagine, as a thought exercise, a large ensemble of macroscopically identical experiments being conducted simultaneously. In this case, *F*_{i} at any time can be interpreted as the instantaneous fraction of the experiments in which the same site *i* inside the macroscopically identical samples remains intact without dislocation generation there. The fraction *F*_{i} should change with time, since as time proceeds, dislocation generation should have occurred at site *i* in more experiments of the ensemble. At any time instant *t*, the nucleation rate at site *i* can be defined as the fractional reduction rate of the number of intact samples within the next time interval d*t*, i.e.(3.2)Equation (3.2) can be integrated with respect to some selected parameter *L* along the (common) strain path in these experiments, which could be the applied force or strain, for example, or even time itself. Such an integration yields(3.3)where . Substituting equation (3.3) into equation (3.1) gives(3.4)In what follows, we assume that the nucleation rate at site *i* scales with the factor *P*(*x*=*x*_{c}) in equation (2.3) for some critical value of *x*_{c}, which can be set to be where *α*≈0.5, i.e. about one-half of the *a*_{o}〈112〉/6 Shockley partial vector. In other words, instability is assumed to occur when the relative displacement of an atomic pair reaches the unstable point for Shockley partial vector formation. An element d*V*_{i} at site *i* in the entire material volume can turn into a nucleus for dislocation generation if any of the atoms inside d*V*_{i} exhibits a relative displacement equal to the critical value, and so the nucleation rate at site *i* is(3.5)where *Ω* is the atomic volume and *ν* is an attempt frequency, which can be taken to be the average atomic vibrational frequency. Substituting equation (3.5) into equation (3.4), we have(3.6)where is the local strain at position inside the specimen at a point in the strain path specified by . In a given loading scheme and material geometry, the strain field is known and hence the probability of survival *F* of the entire material can be evaluated from equation (3.6).

## 4. Application to uniform deformation at constant strain rate

As a special case, the strain field in equation (3.6) can be uniform inside the material volume, as would be the case in uniform compression or shearing. In this case, equation (3.6) becomes(4.1)where *N* is the number of atoms in the material.

### (a) Constant-rate shearing

In a constant-rate shearing experiment, the survival probability of the entire sample can be obtained from equation (4.1) by setting the strain-path variable *L* to be the instantaneous shear strain *γ*,(4.2)Here is the constant strain rate. Equation (4.2) can be used to predict the critical shear strain in the MD simulations described in §2 using the various constants determined from the MD simulations. Figure 5 shows the variation of the survival probability with the instantaneous strain according to equation (4.2), using the following values of the constants determined from the MD simulations at 300 K: *A*=0.174 nm, *B*=0.0182 nm, *m*=3.7, as well as the following constants used in the MD simulations: *N*=31 104, *a*_{o}=0.3573 nm, . Furthermore, *α* is set to be the ideal value of 0.5, and the attempt frequency is assumed to be *ν*=10^{13} s^{−1}, although its precise value is found to have insignificant effects on the results. Figure 5 shows that the survival probability undergoes a rather sharp transition at *γ*≈0.1. For strains below ∼0.08, *F*=1 and so the material is predicted to survive, i.e. be elastic. For strains larger than *γ*≈0.12, *F*=0 which means plastic deformation will definitely happen. It is within a narrow range of *γ*=0.1±0.01 that the material behaviour is predicted to be uncertain, with a finite chance of 0<*F*<1 that the material is elastic. The sharp transitional value of *γ*≈0.1 can therefore be interpreted as the critical applied strain across which the material deformation changes from definitely elastic to definitely plastic. The critical strain observed directly from the load drop or energy relaxation in the MD simulation at 300 K was also 0.1 (see figures 1 and 2), in excellent agreement with the prediction using equation (4.2).

The above calculation was repeated for other temperatures from 20 to 600 K, using the various constants predicted from the MD simulation at each temperature. In figure 6*a*, the theoretical critical strains calculated from equation (4.2), with *α* set to be 0.5, are compared with those observed directly from the MD simulations. It can be seen that at 300 K, the agreement is excellent as discussed above, but the agreement deteriorates as the temperature deviates from this value. However, after slight adjustment in the *α* value at each temperature, the critical strain predicted from equation (4.2) can be brought into agreement with the critical strain observed directly from the MD simulations. Figure 6*b* shows the adjusted values of *α* at different temperatures. It can be seen that the adjusted *α* values fall between 0.4 at 20 K and 0.58 at 600 K, a rather narrow range centred about the expected value of 0.5. The slight deviations of *α* from the ideal value of 0.5 may indicate slight changes in potential barrier along the Shockley partial vector at different temperatures. One may therefore conclude that equation (4.2) can describe fairly satisfactorily the critical strains observed in the MD simulations.

### (b) Uniform compression

Recently, Uchic *et al*. (2004) performed nanoindentation experiments to study the compressive stress–strain responses of 〈123〉-oriented, annealed Ni_{3}Al–Ta micro-pillars within the diameter range of 0.5–20 μm and aspect ratios from 2 : 1 to 4 : 1, and observed a marked increase in the yield point as the specimen size decreases. In a recent paper (Zuo *et al*. 2005), we employed the displacement profile in equation (2.3) to explain the volume dependence of yield strength within a deterministic framework. Here, we re-address the same issue by using the current survival probability consideration. By setting the strain-path variable *L* in equation (4.1) to be the compressive direct stress *σ*, as this is the observable in Uchic *et al*.'s experiments, equation (4.1) can be rewritten as(4.3)Here, *G* is the shear modulus along the incipient slip system, and *M* is a resolution factor connecting this slip system and the compression axis. Hooke's law is assumed in equation (4.3), and this should be approximately valid before the onset of plasticity. For the 〈123〉 orientation used in the experiments by Uchic *et al*. (2004), *M*=1.905 if we assume a slip system. Also, *G*=70.3 GPa for the slip system, calculated using the experimental elastic constants Voter & Chen (1987) used to fit their potentials. The stress rate is not reported by Uchic *et al*. in their paper. Nevertheless, the critical stress predicted from equation (4.3) is found to be very insensitive to the exact value of , and so an indicative value of is assumed here. The number of atoms *N* in a micro-pillar can be calculated as(4.4)where *d* is the micro-pillar diameter, and an aspect ratio of 3 : 1, accounting for the term 3*d*, is assumed here, which is the mean value used in the experiments by Uchic *et al*. The attempt frequency in equation (4.4) can again be assumed to be *ν*=10^{13} s^{−1}. Using the MD simulated values for *A*, *B* and *m* at 300 K, the *F*(*σ*) curve is found to exhibit a sharp transition similar to figure 5 at each micro-pillar size. Figure 7 shows the variation of the critical stress at the sharp transition with the pillar size, using the MD-predicted values of *m*=3.7, *A*=0.174 nm and *B*=0.0182 nm at 300 K, and *α*=0.5. Also shown in figure 7 are the experimental compressive yield strengths reported by Uchic *et al*. It can be seen that, although the theory here predicts the increasing trend of the yield strength as the pillar size deceases, the experimental values are significantly lower than the theoretical prediction. Good agreement between the experiments and the theory, however, can be achieved if the *A* and *B* values are adjusted to 0.086 and 0.02175 nm, respectively. It may be argued that the adjustment is required because the Voter–Chen potential used in the present MD work cannot be expected to very accurately describe the vibrational properties of the ternary Ni_{3}Al–Ta, the material used by Uchic *et al*. (2004), but it is doubtful whether the small ternary addition of Ta can change the potential to such a large extent. It should be noted that the deterministic approach we used earlier (Zuo *et al*. 2005) also significantly overestimated the yield strength if the MD simulated values of *A* and *B* are used. The overestimation is thought to be due to the non-ideal situations in the pillar experiments, which are neglected in the analysis. First, the pillars used by Uchic *et al*. were in the 0.5–20 μm range, and it would be quite unlikely for structures of such sizes to contain no initial defect at all. Secondly, in the pillar geometry, nucleation of dislocations is likely to take place heterogeneously at the pillar free surfaces. Dislocation nucleation from a free surface can be expected to be a lot easier than from the crystal interior. In fact, Horstemeyer *et al*. (2001) used free boundary conditions to simulate the onset of crystal plasticity by MD, and they observed that plasticity was always initiated at the free surfaces of the simulation block. To model the probability of dislocation nucleation from a free surface would require more careful MD simulations and this task is outside the scope of the present paper. However, if nucleation is really from the free surface, the yield condition should also exhibit a similar decreasing dependence upon sample size, because a larger pillar would have more surface atoms and so the chance of nucleation would be higher.

## 5. Application to nanoindentation experiments

### (a) Constant-loading-rate experiments

Recently, the present authors studied incipient plasticity in annealed Ni_{3}Al using nanoindentation (Wo *et al*. 2005). Similar to earlier observations on SiC by Schuh & Lund (2004), in experiments using a constant loading rate, the critical indentation load at which the first strain burst occurred was found to scatter about a modal value in Ni_{3}Al. Here, we attempt to predict this statistical scatter using the present model.

In the indentation situation, the stress field underneath the indenter is not uniform and so the full expression for the survival probability in equation (3.6) has to be used. The strain field underneath the indenter is mathematically very complicated, but since the exponential function in the inner integral in equation (3.6) is rapidly decaying, the inner integral in equation (3.6) is dominated by the maximum value of *γ* at the most highly stressed point, located at a certain depth beneath the indenter along the indentation axis (see p. 94 of Johnson 1985, figure 12). An approximate solution can therefore be obtained by considering only the strain field in the neighbourhood of this most highly stressed point, and to regard this strain field as roughly spherical about the most highly stressed point. Under these approximations, the inner integral in equation (3.6) can be simplified into(5.1)where *r* is a radial coordinate measured from the most highly stressed point, located at a depth of *z*=0.4805 *a* below the surface along the indentation axis, as shown by the analysis in appendix A. From equation (A 5), the strain field *γ*(*r*) in equation (5.1), appropriate for the experimental conditions in our earlier experiments (Wo *et al*. 2005), can be approximately expressed as(5.2)where *a* is the tip–sample contact radius given by equation (A 6), *L* is the indentation load and *G* is the shear modulus. In equation (A 6), the tip radius *R* in our earlier experiments (Wo *et al*. 2005) was determined to be approximately 1.369 μm from direct scanning electron microscopy imaging of the tip, as shown in figure 8. The reduced modulus *E*_{r} is calculated to be 194.2 GPa using the method proposed by Vlassak & Nix (1994) and the elastic constants for Ni_{3}Al commensurate with the Voter–Chen potentials (Voter & Chen 1987), and the known elastic constants for diamond which is the tip material (Young modulus=1140 GPa, Poisson ratio=0.07).

In our earlier experiments (Wo *et al*. 2005), a large number of independent nanoindentation tests were performed on the same grain of the polycrystalline specimen under the same constant-load-rate condition, and the loads at which the strain bursts occurred were recorded to generate a cumulative probability function *F*(*L*) for the proportion of the tests exhibiting the strain burst at a load higher than *L*. Figure 9 shows the probability function *F*(*L*) from the data we reported earlier (Wo *et al*. 2005). Also shown in figure 9 is the theoretical curve predicted by numerically integrating equation (3.6) using equation (5.2) for the inner integral, using the loading rate of =1 mN s^{−1} appropriate for the experimental conditions, as well as the constants where *a*_{o}=0.3573 nm, *ν*=10^{13} s^{−1}, and the MD-predicted constants at 300 K: *m*=3.7, *A*=0.174 nm, *B*=0.0182 nm and *α*=0.509, the adjusted value shown in figure 6*b* at 300 K. The theoretical curve can be seen to be in very good agreement with the experimental results.

### (b) Constant-load experiments

Under a constant load, the nucleation rate in equation (3.5) is constant at a fixed temperature, because the normalization constant *x*_{o} depends only on strain and temperature, according to equation (2.5). From equation (3.2), the survival probability *F*_{i} at each material element should then be an exponentially decaying function with respect to time, with being the decay constant. From equations (3.1) and (3.5), the survival probability of the entire sample can be shown to be(5.3)where the decay constant *κ* is given by(5.3a)The time *t* in equation (5.3) can be regarded as the waiting time before incipient plasticity occurs, and so in a large ensemble of macroscopically identical experiments, equation (5.3) suggests that the probability *F*(*t*) for a given sample to survive longer than time *t* should be exponentially decaying with *t*. The situation is analogous to radioactive decay, for which the nucleation rate is also constant, and the decay time also follows an exponential distribution function.

In our earlier study (Wo *et al*. 2005), we also performed nanoindentation experiments on Ni_{3}Al by holding the load at constant values initially in the elastic regime, and then monitoring the statistical distribution of the waiting time for incipient plasticity to occur. It was found that when the load values were higher than ∼500 μN, the waiting time indeed followed an exponential distribution function. Figure 10 shows a set of data reproduced from our previous paper (Wo *et al*. 2005), from which, it can be seen that the decay constant was about 0.2 s^{−1} at 500 μN. It should be noted that the quantity plotted in figure 10 is the probability density function of the sample lifetime, and this is equal to −d*F*/d*t*, since *F*(*t*) in equation (5.3) is the survival probability, i.e. the probability for a sample to exhibit a lifetime *longer than t*. The theoretical decay constant *κ* in equation (5.3*a*) can be evaluated using the approximate expression for *γ* in equation (5.2). Using the values of *α*=0.509, *m*=3.7, *A*=0.174 nm and *B*=0.0182 nm as used in figure 9, *κ* at 500 μN is calculated to be 245 s^{−1}, which is three orders of magnitude higher than the experimental value of 0.2 s^{−1} in figure 10. However, if *B* is adjusted slightly to 0.017 nm, the correct value of 0.2 s^{−1} for *κ* will be predicted. Alternatively, if the tip radius is not 1.369 μm but is 1.7 μm instead, the correct *κ* value of 0.2 s^{−1} will also be predicted, using the original value of *B*. The precise value of the predicted decay constant is actually very sensitive to small changes in the parameters used in the model, and this is due to the rapidly decaying exponential function in the integrand in the expression for *κ* in equation (5.3*a*). However, the exponential form of the sample lifetime under constant load is correctly predicted using the present model, and the model parameters required to achieve agreement with the experimental decay constant are not inconsistent with those derived from the MD simulations at 300 K.

## 6. Discussion

In the present work, a methodology is proposed to predict the incipient plasticity behaviour of micron-sized crystals from MD simulations on nanometre-sized blocks. The methodology is similar in spirit, but different in nature and details, to the use of Weibull statistics in the prediction of the strengths of large-sized ceramics from laboratory-sized specimens. The Weibull distribution function used in the bulk ceramics context is empirical, but the scalability assumed in the present methodology originates directly from an important assumption in statistical physics—the statistical profiles of atomic properties, such as velocities, energies and so on, are believed to be independent of the system size, as long as it is large enough to enable accurate statistical sampling. In the present methodology, the relative displacement defined in equation (2.2) for each atom is regarded as an atomic property, the statistical profile of which gathered from a nano-sized simulation block is assumed to be the same as that in a larger system. In the entire subject of statistical physics, there is no provision to invalidate such scalability, and indeed, such scalability is a fundamental assumption in statistics. The statistical profile of the relative displacement, *P*(*x*), is in fact a special form of the Boltzmann factor in conventional statistical physics, and so it can be used to predict nucleation rates in much longer time scales than that used in the MD simulation. The *P*(*x*) function in equation (2.3) does not contain any time scale, and its use to predict extremely rare events at strain rates in the experimental regime in fact represents extrapolation far beyond the time scale used in the MD simulations. The present methodology is, therefore, one designed to bridge both time and length scales between the MD and the experimental regimes. However, it should be remembered that the scalability assumed here is valid as long as the material structure does not undergo any systematic or significant change with the system size. It cannot be applied to situations after the generation of the first dislocation, for example, because the presence of dislocations or indeed other crystalline defects will destroy the scalability. Incipient plasticity is perhaps one of the very few problems that can be studied by the present methodology.

Supported by a few experimental findings previously reported in the literature, the present results indicate that micron-sized crystals have the interesting property that their initial yield strength is changeable, or unpredictable in some sense, within a certain load range. Within this load range, one can only speak about a certain probability of survival between 0 and 1, as shown in figure 5 or 9 within the transitional zone. Inside this transitional zone, one knows the *chance* that a component will survive, but not whether a given component will survive or otherwise with absolute certainty. Outside this transitional zone, however, the chance of either surviving or failure will become more certain, as the survival probability transits quickly to 1 or 0, as illustrated in figures 5 and 9. In the design of micro-machinery, this transitional load range should therefore be avoided, or the machine performance will be unpredictable. The range of this transition is also size-dependent from the present analysis. In fact, by differentiating equation (4.2), one finds that(6.1)and so the strain or load range of the transition is approximately proportional to the quantity , and as *N*→∞, the transition range vanishes. Macro-sized specimens are therefore deterministic in terms of their incipient plasticity behaviour, let alone that large specimens should always contain some initial defects which will act as sources for incipient plasticity.

Finally, it can be seen from figures 7 and 9 that, although the general trends of size and load dependence of the survival probability are predicted, the results are very sensitive to the model parameters, such as *m*, *A* and *B*, the values of which can only be obtained by MD simulations. This calls for the use of more reliable atomic potentials in future work in this area. Also, although the present MD simulations are pertinent to only Ni_{3}Al, the general framework of the present methodology should be applicable to other materials, provided that MD simulations are first carried out to identify the incipient plasticity process.

## 7. Conclusions

In this work, MD simulations were carried out to investigate the onset of homogeneous shear deformation in a small block of Ni_{3}Al. Dislocation generation was observed to take place at atomic sites with interatomic displacements approaching one-half of the Shockley partial Burgers vector. The probability distribution function of the interatomic displacements was obtained from the MD simulations, and this was used to derive the survival probability of materials with different volumes supporting different loading conditions. The survival probability was then calculated for several situations including homogeneous deformation and nanoindentation. The results indicated that the survival probability exhibited rather sharp transitions at certain load ranges in all these situations.

## Acknowledgments

This research was supported by a research grant (no. HKU 7201/03E) from the Research Grants Council, Hong Kong Special Administrative Region of People's Republic of China. Helpful discussion with Dr G. P. Zheng is gratefully acknowledged.

## Footnotes

The electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2005.1645 or via http://www.journals.royalsoc.ac.uk.

- Received March 9, 2005.
- Accepted December 14, 2005.

- © 2006 The Royal Society