A Novel Approach to Modelling Nanoindentation Instabilities

Abstract

We review the recently developed models for load fluctuations in the displacement controlled mode and displacement jumps in the load controlled mode of indentation. To do this, we devise a method for calculating plastic contribution to load drops and displacement jumps by setting-up a system of coupled nonlinear time evolution equations for the mobile and forest dislocation densities by including relevant dislocation mechanisms. These equations are then coupled to the equation defining constant displacement rate or load rate. The model for the displacement controlled mode using a spherical indenter predicts all the generic features of nanoindentation such as the elastic branch followed by several force drops of decreasing magnitudes and residual indentation depth after unloading. The stress corresponding to the elastic force maximum is close to the yield stress of an ideal solid. The predicted numbers for all the quantities match experiments on single crystals of Au using a spherical indenter. We extend the approach to model the load controlled nanoindentation experiments that employ a Berkovich indenter. We first identify the dislocation mechanisms contributing to different regions of the $F-z$ curve as a first step for obtaining a good fit to a given experimental $F-z$ curve. This is done by studying the influence of the parameters associated with various dislocation mechanisms on the model $F-z$ curves. The study also demonstrates that the model predicts all the generic features of nanoindentation such as the existence of an initial elastic branch followed by several displacement jumps of decreasing magnitudes and residual plasticity after unloading for a range of model parameter values. Furthermore, an optimized set of parameter values can be easily determined that give a good fit to the experimental load–displacement curves for Al single crystals of $\left(110\right)$ and $\left(133\right)$ orientations. Our model also predicts the indentation size effect in a region where the displacement jumps disappear. The good agreement of the results of the models with experiments supports our view that the present approach can be used as an alternate method to simulations. The approach also provides insights into several open questions.

1. Introduction

The fact that mechanical properties of small volume systems are different from the bulk has been evident in a number of early studies. For instance, the high strength of whiskers is due to the absence of dislocations [1,2]. Similarly, the grain boundary strengthening mechanism is due to the fact that dislocation motion is limited by the grain size, and the dynamic friction (wear) is controlled by the plastic deformation of micrometer or sub-micrometer asperities [3]. There has been a spurt of activity in size dependent studies on plastic deformation of small volume systems in the last three decades due to technological importance as well as the scientific challenges it offers. For instance, intermittent flow is observed when the diameter of micrometer rods are below a certain value while it is smooth when it is large implying the instability manifests when the aspect ratio is reduced [4,5,6]. Similar intermittent plastic flow in the form of load fluctuations or displacement jumps is reported when the indentation depth is less than 100 nm both in thin and bulk samples [7,8,9,10,11,12]. Thus, nanoindentation experiments fall into the class of experiments where plastic instability manifests due to small deformed volume [7,8,9,10,11,12].

Traditionally, two modes of nanoindentation experiments are employed, namely, load/force controlled (LC) mode or displacement controlled (DC) mode [7,9,11,13,14,15]. In both modes, experiments measure the force response of the sample (to the applied force) as a function of indentation depth. While several load drops of decreasing magnitudes are seen in the DC mode experiments [7,9,11,13,14,15], several displacement jumps of decreasing magnitudes are seen beyond the elastic limit in the LC mode experiments [10,12,14,16,17]. The maximum load on the elastic branch is close to the theoretical yield stress. In both modes, residual plasticity is seen after unloading the indenter. While the presence of surface defects, oxide coatings, precipitates, etc. complicate the interpretation of the results, careful and controlled experiments on well-prepared single crystals have demonstrated that the above generic features are reproducible [11,14,15,18]. Since the deformed volume prior to the onset of intermittent plastic flow is a few nm${}^3$, the standard explanation for the high yield strength is the low probability of finding dislocations in such small deformed volumes. The sequence of load drops in the DC mode and displacement jumps in the LC mode beyond the elastic branch have been considered as signatures of instabilities triggered by bursts of plasticity.

Three distinct types of studies have provided much insight into the nanoindentation process, namely, bubble raft indentation[10], colloidal crystals, which are soft matter equivalent of the crystalline phase [19], and in situ transmission electron microscope studies. The latter, in particular, has been useful in visualizing the dislocation nucleation mechanism in real materials [20]. Finally, considerable theoretical understanding has come from various types of simulation studies such as molecular dynamics (MD) simulations [21,22,23,24,25,26], dislocation dynamics simulations and multiscale modeling simulations (using MD together with dislocation dynamics simulations) [26]. While there are a number of MD simulations for the DC mode, there are no simulations of any kind that predict displacement jumps reported in the LC mode experiments.

However, a serious limitation of these simulations is the limited size of the simulated volumes and short time scales used. More importantly, the indentation rates are often several orders of magnitude higher than those in experiments, raising questions of relevance of these simulations to real materials [23,27]. Furthermore, simulation approaches cannot employ experimental parameters such as the indentation rates or the geometrical parameters defining the indenter, thickness of the sample, etc. As a result, the predicted values of load, indentation depth, etc. differ considerably from those reported by experiments.

The above limitations of simulation methods prompted us to develop an alternate framework that has the ability to adopt laboratory time and length scales to predict numbers that agree with experiments. The purpose of this paper is to review the recently introduced dynamical approach and the models that predict the results for the DC and LC mode nanoindentation experiments [28,29]. The approach is based on our previous experience in modeling a significantly more complex spatio-temporal instability, namely, the Portevin–Le Chatelier (PLC) effect seen in bulk samples of dilute alloys under constant strain rate conditions [30,31,32,33,34]. The approach is applicable to ideal single crystals without any kind of defects.

Our approach to nanoindentation is based on two observations. First, the sequence of load drops (in the DC mode) and the sequence of displacement jumps (in the LC mode) seen beyond the elastic branch are triggered by collective dislocation activity. Second, the intermittent deformation is a signature of an instability. The most suitable mathematical platform for describing instabilities and collective behavior of constituent entities is nonlinear dynamics. Furthermore, since the basic defects contributing to plastic flow are dislocations, a proper description of collective dislocation behavior demands that we include all the relevant dislocation mechanisms contributing to the bursts of plasticity. On the basis of these observations, we recently developed a novel approach [28,29] by combining the power of nonlinear dynamics with a method for calculating the plastic contribution to the indentation depth. The latter is calculated by developing time evolution equations for the mobile and forest dislocation densities based on the knowledge of the dislocation mechanisms contributing to plastic flow [28,29]. In addition, the approach allows us to employ experimental indentation rates, the geometrical parameters defining the indenter, thickness of the film, etc. The efficacy of the approach is illustrated by developing a model for the DC mode nanoindentation experiments (carried out using a spherical indenter) [28]. The model predicts all the generic features of the DC mode with numbers matching experimental results on single crystals of gold [11,15]. We have also developed a model for predicting the displacement jumps in the LC mode indentation even though the exercise is even more challenging due to the conceptual difficulty in enforcing a constant load rate during displacement jumps. Indeed, this difficulty appears to be the reason for a complete absence of any kind of simulations to explain the displacement jumps in the LC mode. Since we plan to get a good fit to the experimental $F-z$ plots [12], we first study the influence of the parameters associated with dislocation mechanisms on the model $F-z$ curves with a view to identify the dislocation mechanisms that control different regions of the $F-z$ curve. The study further demonstrates that the model predicts all the generic features of nanoindentation such as the existence of an initial elastic branch followed by several displacement jumps of decreasing magnitudes and residual plasticity after unloading for a range of values of model parameters associated with the dislocation mechanisms. This also helps us to obtain optimized parameter values that give a good fit to the experimental load–displacement curves for Al single crystals of $\left(110\right)$ and $\left(133\right)$ orientations [12]. While hardness is not well defined at nanometer scales where displacement jumps dominate, for larger indentation depths where intermittency disappears, the hardness decreases with depth. Thus, indentation size effect is also predicted. These results also justify our expectation that the present approach can be used as an alternate method to simulations in modeling nanoindentation instabilities. These studies also provide insights into several open questions related to the collective behavior of dislocations: (a) Which dislocation mechanisms determine the first and subsequent load drops or displacement jumps? (b) Can all load drops and displacement jumps be explained on the basis of a common physical and mathematical mechanism? (c) What physical and mathematical mechanisms contribute to the decreasing magnitudes of the load drops or displacement jumps?

2. Insights from Simulations on Nanoindentation

At the outset, it must be mentioned that, apart from bubble raft simulation [10], only MD simulations predict load drops in DC mode of nanoindentation. There are no MD simulations (or any other for that matter) for the displacement jumps in the LC mode. In general, simulation methods have always served the purpose of understanding dislocation mechanisms underlying the evolution of complex dislocation microstructures in various dislocation mediated plasticity. The purpose of this section is to briefly summarize the results of MD simulations and multiscale modeling method, and discuss their advantages and limitations.

In the context of simulating the DC mode nanoindentation instability, MD simulations have played an important role in providing a good understanding of initial stages of plastic deformation. However, because of the limited volume that can be simulated, subsequent evolution of the plastic zone becomes prohibitive as indentation proceeds. This has motivated researchers to develop a multi-scale modeling method that uses MD results as a starting point for understanding further evolution of dislocation microstructure using discrete dislocation dynamical (DDD) simulations in conjunction with finite element methods (FEM).

The suitability of these simulations depends on the stage of evolution of dislocation microstructure that is being targeted. Here, we note that the load drops manifest at a few tens of nm scale corresponding to the initial stages of dislocation microstructure. Furthermore, the initial indented volume prior to the onset of plasticity is of the order of 10 nm${}^3$. On the other hand, typical dislocation density in a well grown single crystals is ∼10${}^{10}/$m${}^2$. Thus, the probability of finding a dislocation is close to zero. The nucleated dislocation loops act as the trigger for the pop-in events and therefore represent elastic to plastic transition. Then, MD simulations are best suited for this initial stages of deformation. The typical size of the simulated volume is a few million atoms. However, they are clearly not suitable for simulations when describing further development of plastic zone (beyond say 50 nm). A multi-scale method has been developed by Chang et al. [26] that incorporates MD results into discrete dislocation dynamic(DDD) simulations in conjunction with finite element methods to take into account the boundary conditions relevant for the loading condition.

There are a number of MD simulations that capture the evolution of dislocation microstructure to various degrees [21,23,24,25,26]. Features that emerge from these simulations include the nucleation of dislocation loops, their expansion and interaction, detachment from the source, and their propagation [26].

Almost all studies use a spherical indenter, which is modeled using a repulsive potential. Different kinds of potentials have been used (to simulate the desired crystal) in MD simulations. Typical large MD simulations use a million atoms and the radius of the indenter can range from about 0.8 nm to 15 nm [21,23,24,25,26]. Typical rates of indentation are ∼$0.1\mathsf{\AA }$/step with each step typically a few ps duration required to attain minimum energy configuration. Typical indentation depths are about 10 Å and the force is ∼$0.1\mathsf{\mu }$N to a few $\mathsf{\mu }$N.

For illustration, we use the results of Van Vliet et al. [23]. These authors use Ercolessi–Adamssi potential to simulate Aluminum. The load–displacement curve with several load drops is shown in Figure 1a. Due to the inherent advantage of MD simulations, dislocation activity (see Figure 9 of Ref. [23]) can be correlated with the load drops in the load–displacement curve. As can be seen, nucleation occurs at 20 nN, which corresponds to a nucleation stress of ${\sigma }_m$∼ 3.4 GPa. The glide loops expand along $\left\{111\right\}$ planes. In addition to nucleation, other dislocation mechanisms the authors realize are the formation of sessile locks, cross-slip, and propagation after detachment from the source. The maximum depth of indentation is $20\mathsf{\AA }$. One can also notice small scale fluctuations in the load–displacement curve, an effect well known in MD simulations, even in the context of MD simulations of equilibrium properties. Such fluctuations tend to smooth out as the system size is increased.

In the following, we briefly summarize the results of the multi-scale modeling approach due to Change et al. [26]. To begin with, we briefly recall their MD results. They use an Embedded-atom method type of interatomic potential typical to several other MD simulations. The authors use a repulsive sphere that is pressed into the sample at a rate of $0.1\mathsf{\AA }$ at each step of a few ps duration. The top surface is traction-free and the atoms located on the bottom surface are fixed. Periodic boundary conditions are applied in the direction perpendicular to the indentation axis. Figure 1b shows the load–displacement curve with several pop-in events. The first load drop has been identified with the nucleation of three prismatic interstitial dislocation loops. The nucleation occurs at 1 nm and subsequent evolution of dislocation microstructure with indentation depth consists of expansion of the loops, their multiplication followed by detachment from the surface (see Figure 2a–f of Ref. [26]).

The DDD simulations employed by the authors use an Edge-screw code with an FEM solver. The FEM is used to compute the highly heterogeneous stresses during the indentation process. The information about nucleation of dislocation loops seen in MD simulation are incorporated by a set of rules that specify the shape of loops and position of the loops to be introduced in the simulation cell. The imposed displacement rate determines the number of loops introduced. The load–displacement curve using a spherical indenter of radius $R=150$ nm is shown in Figure 2. The figure shows the total load, load contribution from dislocations and the MD result extrapolated to the maximum depth of 100 nm. The evolution of dislocation microstructure becomes increasingly complex as illustrated in Figure 3a–d. For details, see Chang et al. [26].

As can be seen, a major advantage of the MD simulations is their ability to include a range of dislocation mechanisms starting from the nucleation of dislocations, their multiplication, formation of locks, junctions, propagation of loops, etc. [21,23,24,25,26]. Note that these dislocation mechanisms are all used in our model. The most obvious limitation of MD simulations is the limited volume they can simulate. This also means that the maximum radius of a spherical indenter is also limited. This also means that the predicted depths at which nucleation occurs is typically 6–10 Å significantly smaller than $65\mathsf{\AA }$ in experiments [11,15]. One other important limitation is the length and time scales that are inherent to the small scale simulations. The relaxation time scale to attain the minimum energy configuration is typically a few ps. This coupled with the fact that the loading rates are a fraction of an Å means that the imposed displacement ∼m/s. This is at least ${10}^8$ orders higher than that used in experiments by Kiely et al., who use $0.5\mathsf{\AA }/\mathrm{s}$ [11,15]. On the other hand, while DDD simulations are useful in tracking further development of dislocation microstructure, these simulations do not model the intermittent region. These limitations prompted us to devise an alternate framework that can use experimental conditions such as the imposed rates, geometrical parameters defining the indenter and other experimental parameters.

3. Background Material

This review specifically targets two most popular quasi-static modes of indentation, namely, the displacement controlled indentation where imposed displacement is increased at a constant rate and load controlled mode where the load is increased at a constant rate. These two modes of indentation are primarily chosen because the equations governing the DC and LC modes take a simple form compared to the dynamic loading mode, not considered here. The simplicity of these equations helps us establish a good correlation between imposed rates and plastic response of the material. In addition, for each of these modes of indentation, we can use different types of indenters. Here, we consider (a) spherical indenter for the DC mode with a view to compare the model results with those of Kiely et al. [11,15]), and (b) Berkovich indenter in the LC mode so as to compare our model results with Gouldstone et al. [12]. Here, we consider a single crystal without surface defects, inclusions, precipitates, etc. (Other complications in nanoindentation such as pile-up and sink-in effects are beyond the scope of our approach.). We begin by collecting relevant background material.

In our approach, the indentation depth z measured from the undeformed surface $z=0$ is a dynamical variable, i.e., a variable that evolves with time. The total indentation depth z is the sum of elastic indentation depth $z_e$ and plastic depth $z_p$. Then,

$$z=z_e+z_p.$$

(1)

The elastic depth is essentially known since it is governed by the indenter specific load–depth expression. However, calculating the plastic displacement $z_p$ has been the obstacle for developing any nanoindentation model. The strength of our approach is that it provides a way of calculating $z_p$ by first calculating the plastic strain rate ${\dot{ϵ}}_p\left(t\right)$ using the Orowan equation [35]. This also means that we use strain as a fundamental variable that can be calculated. Noting that plastic deformation occurs within the sample of thickness T, we define the strain variable $ϵ=\frac{z}{T}$. Then, we have the total strain $ϵ={ϵ}_e+{ϵ}_p$, where ${ϵ}_e$ and ${ϵ}_p$ refer, respectively, to the elastic and plastic components of the strain. Taking the time derivative, we have

$$\dot{ϵ}=\frac{\dot{z}}{T}=\frac{{\dot{z}}_e}{T}+\frac{{\dot{z}}_p}{T}={\dot{ϵ}}_e+{\dot{ϵ}}_p.$$

(2)

Then, Equation (1) takes the form

$$z=z_e+T\int {\dot{ϵ}}_P\left(t\right)dt.$$

(3)

Note that z is well defined for any T. The plastic strain rate itself is given by the Orowan expression [35]

$${\dot{ϵ}}_p\left(t\right)=bV\left(\sigma \right){\rho }_m\left(t\right).$$

(4)

Here, b is the magnitude of the Burgers vector, $V\left(\sigma \right)$ is the mean velocity of dislocations, ${\rho }_m$ the mobile dislocation density and $\sigma$ the stress. It is clear that this requires calculating ${\rho }_m$ as a function of time. This is obtained by developing the time evolution equations for the mobile ${\rho }_m$ and forest ${\rho }_f$ densities and coupling them to equations defining the DC or the LC modes of indentation.

Note that Equation (4) is a function of stress. Therefore, $\sigma$ must be expressed as a function of the load F and area A supporting the load since both evolve with time during indentation. This requires expressions for the load $F\left(z\right)$ and area $A\left(z\right)$ as a function of z. Both are indenter specific.

3.1. Displacement Controlled Indentation Using a Spherical Indenter

Since we plan to address nanoindentation experiments on single crystals of Au carried out in the DC mode using a spherical indenter [11,15], we first recall a few equations. Consider a spherical indenter of radius R. Then, force response of a sample is given by the Hertzian expression

$$F=\frac{4}{3}{E}^{\ast }{R}^{1/2}{z}_e^{3/2}.$$

(5)

Here, $E^{\ast }$ is the effective modulus of the indenter and the sample given by

$$\frac{1}{E^{\ast }}=\frac{1-{\nu }_s^2}{E_s}+\frac{1-{\nu }_i^2}{E_i},$$

(6)

where $\nu$ and E refer to the Poisson’s ratio and Young’s modulus of the sample s and indenter i.

Pile-up and sink-in effects are routinely seen in indentation experiments. In terms of material property, pile-up or sink-in effects depend on the ratio of the elastic modulus to the yield tress and strain hardening property of the material [13]. Pile-up is characterized by larger contact area compared to the ideal case while a sink-in effect is characterized by a smaller contact area. Thus, $\frac{z_e}{z_{e_c}}$, the ratio of the indentation depth $z_e$ measured from $z=0$ to the contact depth $z_{e_c}$, is taken as a measure of pile-up or sink-in effects. Thus, area–depth relations are usually expressed in terms of contact depth $z_{e_c}$ measured from the contact point.

The area expression for a spherical indenter is given by

$$A_H\left(z_e\right)=\pi R{z}_{e_c}.$$

(7)

Note that Equations (5) and (7) are equations obtained in mechanical equilibrium and therefore do not have any time dependence at this point. However, once these equations are used in equations governing the loading condition and evolution equations for dislocation densities, they acquire time dependence.

Now, consider a spherical indenter of radius R driven into a sample of thickness T at a rate $\dot{r}$ (m/s). Rearranging Equation (2) and substituting, $\dot{z}=\dot{r}$ and ${\dot{z}}_p=T{\dot{ϵ}}_p$, the machine equation corresponding to the DC mode nanoindentation reads

$${\dot{z}}_e=\dot{r}-T{\dot{ϵ}}_p.$$

(8)

Note that this equation is applicable to any indenter. For a spherical indenter of radius R used by Kiely et al. [11,15], we can express $z_e$ in terms of F. Then, Equation (8) takes the form

$$\frac{d}{dt}{[\frac{3F}{4E^{\ast }{R}^{1/2}}]}^{2/3}=\dot{r}-T{\dot{ϵ}}_p.$$

(9)

This equation is the force response of a sample subjected to constant displacement rate.

3.2. Load Controlled Indentation Using a Berkovich Indenter

Now, consider the LC mode nanoindentation using a Berkovich indenter with a view to compare our model results with the LC mode experiments on single crystals of Al that employs a Berkovich indenter [12]. Since our equations are functions of stress, the first step is to adopt phenomenological expressions for the load the area supporting the load as functions of indentation depth for a Berkovich indenter. Here, we use the phenomenological expression for the load [12] given by

$$F=2.189E^{\ast }[1-0.21{\nu }_s-0.01{\nu }_s^2-0.41{\nu }_s^3]z_e^2=C{E}^{\ast }{z}_e^2.$$

(10)

Unlike Equation (7), the area function for a non-ideal Berkovich indenter is complicated due to the complicated geometry of the blunted the tip of the indenter [36,37,38,39]. Several expressions have been suggested in the literature for the area of a non-ideal Berkovich indenter. The most general form is the power law expression $A={\sum }_{n=0}{c}_n{z}^{2/2^n}$ [37]. The calibrated shape of the tip of the indenter can be approximated by a spherical shape with a nominal radius R determined by SEM measurements [23,38,39]. Here, we follow Bei et al. The area of a blunted Berkovich indenter is taken to be sum of the area of an ideal Berkovich indenter and the area of a spherical indenter of radius R. Using the area of a sharp Berkovich indenter given by $A_B=24.54z_{e_c}^2$, the area function of a real Berkovich intender is given by

$$A\left(z_{e_c}\right)=24.54z_{e_c}^2+2\pi R{z}_{e_c}.$$

(11)

Bei et al. [38] show that Equation (11) fits the measured area closely for $C{r}_{3}Si$ single crystals. The geometry of the indenter also provides a relationship between $z_e$ and the depth from the contact point $z_{e_c}$, i.e., $z_e=z_e\left(z_{e_c}\right)$ [38]. A linear relation between $z_e$ and $z_{e_c}$, i.e., $z_e=s{z}_{e_c}$ holds for the range of indentation depths used in nanoindentation studies.

While Equations (7) and (11) strictly hold in the elastic region only, we assume that they hold over the entire duration of indentation that includes the elasto-plastic region. i.e., $A\left(z_e\right)\to A(z_e+z_p)$. Indeed, in MD simulations, the measured area increases in steps at each burst of plasticity (see [23]). This assumption is further supported by the fact that the load increases after every load drop or displacement jump. This can only happen if the area supporting the load increases abruptly to make contact with the indenter surface after every load drop or displacement jump.

For the LC mode of indentation, the equation for a constant load rate is given by

$$\frac{dF\left(z\right)}{dt}={\dot{F}}_0=\mathrm{constant},$$

(12)

where ${\dot{F}}_0$ is the constant load rate. (In experiments, the maximum desired load is achieved in a predetermined number of steps.) This equation is valid for any indenter. For the Berkovich indenter, the expression for the load given by Equation (10) is used.

It is worthwhile to comment on the possible influence of the pile-up or sink-in effects on the properties addressed in present context of intermittent indentation process. It is clear that, since the area supporting the load is larger (pile-up) or smaller (sink-in) than the ideal case and the area expression is a function of $z_{e_c}$, the pile-up or sink-in effects are incorporated in the scale factor, s in $z_e=s{z}_{e_c}$. However, the measured area (or equivalently the scale factor s) is seldom given in experiments on nanoindentation instabilities. Thus, one needs to use a judicious choice of s. Since the focus of the work is in formulating a theoretical basis for describing the two (DC and LC mode) nanoindentation instabilities, any judicious choice of the scale factor, s that accounts for all aspects relating to intermittent indentation process is adequate. We shall later comment on hardness calculation in the two modes of nanoindentation instabilities. However, since pile-up or sink-in effects require spatial degrees of freedom for their characterization, calculating pile-up or sink-in effects independently is beyond the scope of our approach since our approach does not include spatial degrees of freedom.

4. Dynamical Approach to Modeling Nanoindentation Instabilities

The basic premise of our approach is that specimen or volume averaged dislocation densities are adequate to describe the measured load-indentation depth curves. This is based on the observation that experimental load–displacement curves are also specimen averages of the dislocation activity in the sample. For the same reason, experimental $F-z$ curves do not have any information about the inhomogeneous nature of the plastic deformation occurring within the sample. Thus, as long as we include all relevant dislocation mechanisms contributing to bursts of plasticity, it should be possible to model the time development of the load and indentation depth.

Further support for this premise comes from the fact that similar time evolution equations for the volume averaged dislocation densities have been effectively used in the literature to explain the temporal aspects of several spatio-temporal plastic deformation instabilities and patterns [30,32,34,40,41,42,43]. For example, the ’storage recovery’ model for the evolution of the forest dislocation density (with respect to shear strain or time) has been used to obtain good insight into several characteristic features of dislocation cell pattern seen in stage III deformation of F.C.C materials [40]. The equation successfully predicts a number of experimental results such as the saturation stress, the hardening rate in the stage II, etc. Similarly, the temporal aspects of the PLC effect manifesting in the form of three types of serrations found with increasing strain rate has been well captured by the Ananthakrishna (AK) model that uses three types of volume averaged dislocation densities (see, for instance, Refs. [30,32,41,42]). Similar volume averaged dislocation densities have also been used to explain the stress–strain curve for a smooth yield phenomenon in the context of acoustic emission [34]. The same type of approach has been used to explain the inverse power law dependence of the yield stress on the diameter of nano-pillars [43]. All these phenomena are spatio-temporal in nature, yet using time dependent equations have proved quite effective in predicting the salient features. While this kind of mean field approach is not common in the plasticity area, it is routinely used in several areas of physics, chemistry and engineering [44,45].

General Form of Time Evolution Equations for Dislocation Densities

Following our earlier publications [28,29,30,32,33,34,41,46], we develop the time evolution equations for the mobile ${\rho }_m$ and forest ${\rho }_f$ densities based on relevant dislocation mechanisms contributing to the plastic flow. Dislocation mechanisms can be broadly categorized into dislocation production and transformation mechanisms. The former includes dislocation nucleation and multiplication. (In our model, all mechanisms leading to the line length increase are regarded as dislocation multiplication processes). Furthermore, since the sample is expected to be dislocation free within the initial indented nanometer volume, nucleation of dislocation loops has to occur before their multiplication. A few simulations [16,21,23,26] also show detachment of the loops from the source followed by their propagation to the boundary (see Figure 9 of Ref. [23] and Figure 5 of Ref. [26].). Each of these dislocation mechanisms can only be activated beyond a certain threshold stress. We denote the threshold stresses for nucleation, multiplication and propagation, respectively, by ${\sigma }_n,{\sigma }_m$ and ${\sigma }_p$. The rate of nucleation of dislocations per unit area is given by $\frac{{\nu }_n}{\pi b^2}exp(\sigma /{\sigma }_n)$. For the sake of simplicity, we assume that a single loop is nucleated, i.e., $\frac{{\nu }_n}{\pi b^2}=1$.

The rate of multiplication of dislocations is traditionally written as $\theta V_m\left(\sigma \right){\rho }_m$, where $V_m\left(\sigma \right)$ represents the average velocity of dislocations [47] and $\theta$ the inverse of an appropriate length scale [28,30]. Commonly used phenomenological expressions [28,29,48] are

$$V_m\left(\sigma \right)=V_{0}exp[\frac{\sigma -h{\rho }_f^{1/2}}{{\sigma }_m}],$$

(13)

and

$$V_m\left(\sigma \right)=V_0[\frac{\sigma -h{\rho }_f^{1/2}}{{\sigma }_m}]{}^m.$$

(14)

In Equations (13) and (14), $h{\rho }_f^{1/2}$ is the back stress with $h=\alpha Gb$, $\alpha$∼0.3 and G is the shear modulus. In Equation (14), m is the velocity exponent. We shall use Equation (13), for the DC mode experiments using a Spherical indenter and Equation (14) for the LC mode where a Berkovich indenter is used.

In time sequence, detachment of dislocations from their source and their subsequent propagation to the boundary follows dislocation multiplication. The rate of propagation of dislocation loops is proportional to the time required for dislocations of velocity $v_p\left(\sigma \right)$ to reach the boundary taken to be located at $z=T$. Then, ${\nu }_p=\frac{v_p}{T}$. We shall also assume that $v_p\left(\sigma \right)$ obeys a power law of the form $v_p=v_{p_0}{[\frac{\sigma }{{\sigma }_p}]}^p$, where p is an exponent. Then, ${\nu }_p=\frac{v_{p_0}}{T}{[\frac{\sigma }{{\sigma }_p}]}^p={\nu }_{p_0}{[\frac{\sigma }{{\sigma }_p}]}^p$. Note that $v_p\left(\sigma \right)$ is the propagation velocity of dislocations through an undeformed part of the specimen. Therefore, p is necessarily different from the exponent m in $V_m\left(\sigma \right)$.

As dislocations multiply, they begin to interact with each other to form junctions. As discussed in a number of earlier publications, two distinct types of short range interactions can be identified that act as the source for the growth of the forest density ${\rho }_f$. These act as loss terms to the mobile density ${\rho }_m$ [28,30,33,34,42,49]. When two dislocations moving in nearby glide planes approach a minimum distance (typically a few nanometers), they can form dipoles. The corresponding loss term to ${\rho }_m$ is $\beta {\rho }_m^2$. Then, $\beta {\rho }_m^2$ is a source (storage) term for the growth of ${\rho }_f$. Here, $\beta$ is a ’rate constant’ that has the dimension of the area covered per unit time. A mobile dislocation can annihilate a forest dislocation. This is represented by $f\beta {\rho }_m{\rho }_f$. The parameter f is a dimensionless constant typically ∼${10}^{-2}$–10. This is a common loss term for both ${\rho }_m$ and ${\rho }_f$. Finally, unlike the dipole source term that has a fixed short range interaction of a few nanometers, there is another source term for ${\rho }_f$ whose interaction range evolves with deformation. This happens when dislocations moving in different glide planes intersect each other to form junctions. This can be written as $\mathsf{\Lambda }{\rho }_m{\rho }_f$. Here, $\mathsf{\Lambda }$ refers to the mean separation between the junctions that also evolves with deformation (or time) since ${\rho }_f$ itself evolves with time. Noting that the mean distance between the forest dislocations is proportional to $1/{\rho }_f^{1/2}$, we may write $\mathsf{\Lambda }=\delta {\rho }_f^{-1/2}$. Here, $\delta$ is taken to be a constant. Then, the loss term for ${\rho }_m$ is $\delta {\rho }_m{\rho }_f^{1/2}$. Clearly, this is a gain term to ${\rho }_f$. ($\delta$ has the dimension of velocity.) This storage term represents the forest mechanism. The storage and recovery terms ($\delta {\rho }_m{\rho }_f^{1/2}$ and $f\beta {\rho }_m{\rho }_f$, respectively) are the only mechanisms used in the storage-recovery model (except for the difference in the notation) [40]. Collecting the loss and gain terms for ${\rho }_m$ and ${\rho }_f$, the time evolution equations can be written as

$$\begin{array}{ccc}\hfill {\dot{\rho }}_m& =& \frac{{\nu }_n}{\pi b^2}exp\frac{\sigma }{{\sigma }_n}+\theta V_m\left(\sigma \right){\rho }_m-{\nu }_{p_0}{\rho }_m{\left[\frac{\sigma }{{\sigma }_p}\right]}^p-\beta {\rho }_m^2-\delta {\rho }_m{\rho }_f^{1/2}-f\beta {\rho }_m{\rho }_f,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill {\dot{\rho }}_f& =& \beta {\rho }_m^2-f\beta {\rho }_m{\rho }_f+\delta {\rho }_m{\rho }_f^{1/2}.\hfill \end{array}$$

(16)

The first two terms in Equation (15) refer respectively to the nucleation and multiplication of dislocations. The third term represents loss to the boundary due to propagation of dislocations. The next three terms are the dislocation transformation terms. Note also that, unlike the common loss term $f\beta {\rho }_m{\rho }_f$ in Equation (15) and (16), $\beta {\rho }_m^2$ and $\delta {\rho }_m{\rho }_f^{1/2}$ are the storage terms for the forest density that contribute to the growth of ${\rho }_f$. Thus, a competition between the parameters $f\beta$ corresponding to the common loss term, and the parameters corresponding to storage terms, $\beta$ and $\delta$, control the relative magnitudes of ${\rho }_m$ and ${\rho }_f$.

5. A Dislocation Dynamical Model for Load Fluctuations Using a Spherical Indenter

We now recast the general form of Equations (15) and (16) suitable for the DC mode indentation using a spherical indenter [28] and compare the predicted results with the results for single crystals of gold [11,15]. This means that we use Equation (13) for $V_m\left(\sigma \right)$ and express $\sigma =F/A_H$ in Equations (15) and (16) using Equations (5) and (7). A few comments are in order at this point. Indentation is carried out on Au samples that have high $E^{\ast }/{\sigma }_y$ ratio. Therefore, one expects pile-up effects are seen. However, micrograph Figure 1b of Ref. [11,15] shows a mild pile-up effect. The pile-up effect is reflected in the increased area (with respect to an ideal spherical indenter) given by $A_H=\pi R{z}_{e_c}$ in terms of contact depth $z_{e_c}$. However, no information is available about the measured area or the scale factor s in $z_e=s{z}_{e_c}$. Moreover, intermittent indentation process is seen until the maximum depth of 25 nm where hardness is ill defined as we shall see. For these reasons, we have used $s=1$ in the following calculations. However, it must be stated that the basic results and conclusion remain valid for any s.

Then, the evolution equations for ${\rho }_m$ and ${\rho }_f$ for the Hertzian indenter of radius R take the form

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{\rho }_m}{dt}}& \hfill =& \frac{{\nu }_n}{\pi b^2}exp\frac{\frac{F}{\pi R[z_e+z_p]}}{{\sigma }_n}+\theta V_0{\rho }_{m}exp\frac{\frac{F}{\pi R[z_e+z_p]}-h{\rho }_f^{1/2}}{{\sigma }_m}-\beta {\rho }_m^2-f\beta {\rho }_m{\rho }_f-\delta {\rho }_m{\rho }_f^{1/2}\hfill \\ & \hfill -& {\nu }_{p0}{\rho }_m{\left[\frac{\frac{F}{\pi R[z_e+z_p]}}{{\sigma }_p}\right]}^p,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\frac{d{\rho }_f}{dt}\hfill & =\hfill & \beta {\rho }_m^2-f\beta {\rho }_m{\rho }_f+\delta {\rho }_m{\rho }_f^{1/2}.\hfill \end{array}$$

(18)

These equations are coupled to the machine equation (Equation (9)) expressed in terms of load and area. Then, we have

$$\frac{d}{dt}{\left[\frac{3F}{4E^{\ast }{R}^{1/2}{T}^{3/2}}\right]}^{2/3}=\frac{\dot{r}}{T}-b{\rho }_m{V}_{0}exp\{\frac{F}{\pi R[z_e+z_p]{\sigma }_m}\}.$$

(19)

Equations (17) and (18) and Equation (19) constitute a coupled set of nonlinear equations for the nanoindentation problem.

5.1. Estimation of Parameters

Since there are several parameters, we begin by summarizing the estimated values of all the parameters detailed in Refs. [28,29]. Note that the model parameters can be classified as experimental and theoretical. As mentioned in the Introduction, since we work at laboratory (length and time) scales, we can directly adopt experimental parameters such as indentation rate $\dot{r},R,T,E^{\ast },b$ and h. For the present case, we use the experimental parameters used in Ref. [11] corresponding to the gold samples. For this case, the imposed experimental displacement rate is $0.5\mathsf{\AA }/\mathrm{s}$ while the imposed rate in MD simulations are several orders higher [21,22,23,24,25,26].

The model has several theoretical parameters. A comment is desirable in this context. In dislocation based plasticity models, the number of parameters is determined by the number of dislocation mechanisms required to model the phenomenon. Thus, the more complex the phenomenon, the more the number of dislocation mechanisms and the more the number of parameters. For instance, the AK model [30,32,33,46,50] for the complex spatio-temporal PLC instability has ten parameters. Even in simple models such as the storage recovery model [40] and the model for plastic deformation of micro-pillars [43], there are three and six respectively. In the present case, we have ten parameters corresponding to six dislocation mechanisms [28,29].

Now, consider estimating the theoretical parameters $\theta V_0,f\beta ,\delta$ and ${\nu }_{p_0}$. Estimation of the parameter values for the two nanoindentation models has been carried out in detail in Refs. [28,29]. Very briefly, $\theta V_0$ is the basic time scale for the growth of the mobile dislocation density that determines the plastic flow. We set $\theta V_0=1$ to ensure that it matches the experimental time scale. (Here we use $V_0$∼${10}^{-6}$ m/s and $\theta$∼${10}^6/$m). The three threshold stresses ${\sigma }_n,{\sigma }_m$ and ${\sigma }_p$ are fixed based on the basis of other studies or on physical grounds. The nucleation threshold is known to be ${\sigma }_n$∼$E^{\ast }/10$. Similarly, the threshold for multiplication of dislocations is typically ${\sigma }_m$∼$E^{\ast }/300$. As for the propagation threshold ${\sigma }_p$, since the propagation can only occur after dislocation multiplication, we assume ${\sigma }_m<{\sigma }_p$. The upper limit of the velocity prefactor $v_{p_0}$ in the expression for propagation velocity is limited by $V_0$. While the exponent value p is not known, we fix it to be $p=1$.

As for the parameters $f\beta$ and $\delta$, we have demonstrated that the orders of magnitudes of these two parameters are determined once the asymptotic (long time) values of the two densities are given or vice versa [28,29]. See Ref. [28,29] for details.

A further requirement on the allowed ranges of parameter values is that Equations (17)–(19) should exhibit an instability. Since these equations form a coupled set of nonlinear equations, such equations often exhibit instabilities for a domain of values of the model parameters. As demonstrated in Ref. [28], the stability matrix corresponding to the steady state of Equations (17)–(19) exhibits a pair of complex conjugate roots with a negative real part. Thus, the DC mode nanoindentation is a transient instability rather than a true instability. Indeed, this is the underlying mathematical mechanism for the decreasing magnitudes of the load drops. (See also Appendix of Ref. [28]). The instability domain of parameter values has been determined numerically. The analysis shows that the instability occurs over a wide range of values of the parameters [28]. In general, the physically relevant values of the parameters form only a small subset of the instability domain. The precise values of the parameters used in the our calculation are given in the

Table 1. Parameter values used for the model. The experimental parameters used are drawn from experiments on single crystals of Au [11].

${\mathit{E}}^{\ast }$(GPa)	${\mathit{\sigma }}_{\mathit{n}}$(GPa)	${\mathit{\sigma }}_{\mathit{p}}$(GPa)	${\mathit{\sigma }}_{\mathit{m}}$(GPa)	${\mathit{v}}_0$(m/s)
75	$7.5$	$1.91$	$0.955$	${10}^{-6}$
$\dot{\mathit{r}}$ (Å/s)	$\mathit{R}=\mathit{T}$ (Å)	$\mathit{\delta }$ (m/s)	$\mathit{\beta }$ (m${}^{\mathbf{2}}$/s)	$\mathit{f}$
$0.5$	1075	$3.92\times {10}^{-6}$	${10}^{-13}$	$2.5$

.

5.2. Results

Equations (17)–(19) have been solved using an adaptive time step Runge–Kutta solver (ODE15S MATLAB, (MathWorks, Natick, MA, USA)) with initial conditions ${\rho }_m=0$ and ${\rho }_f=0$ consistent with the absence of dislocations in the indentation volume. Using the parameters given in Table 1, we calculate the load as a function of indentation depth. A plot of $F-z$ is shown in Figure 4a. The loading run is shown by the blue curve and the unloading run by the red curve. It is clear that a large load drop is seen after the initial Hertzian response. The load maximum on the Hertzian branch is $F=15.1$ $\mathsf{\mu }$N and the first load drop is ∼10.6 $\mathsf{\mu }$N. These values are close to the experimental values of 15 $\mathsf{\mu }$N and 10 $\mathsf{\mu }$N, respectively [11]. Subsequent load drops show a decreasing trend as in experiments. The corresponding stress as a function of z is shown in Figure 4b. The maximum stress on the Hertzian branch is ≈7.5 GPa (see Figure 4b), the same as used in experiment [11].

Our approach automatically allows us to calculate ${\rho }_m$ and ${\rho }_f$ as a function of time or as a function of depth. Two features are evident from the plot of ${\rho }_m-z$ shown in Figure 5a. First, the nature of ${\rho }_m$ is of burst type with the bursts occurring at each load drop with successive bursts decreasing rapidly from an initial high value of ∼7 × 10${}^{13}$/m${}^2$. Second, the duration of the bursts is very short with the durations of the successive bursts increasing concomitantly. The large magnitudes of ${\rho }_m$ bursts together with their short durations are the characteristic features of collective behavior of dislocations. Another conclusion that can be drawn is that collective effects become weaker with successive load drops and are eventually lost as is clear from the decreasing magnitudes and increasing durations of the bursts of ${\rho }_m$. In contrast, the corresponding forest density ${\rho }_f$ increases in steps (not shown, see Figure 2b of [28]). The magnitude of the steps decrease with z, eventually reaching a near saturation level. The positions of these steps are found to be well correlated with the positions of the ${\rho }_m$ bursts. The stepped increasing nature of ${\rho }_f$ and their decreasing magnitudes is due to the fact that a large part of ${\rho }_m$ bursts are transformed to ${\rho }_f$. This also implies depletion of the peak heights of ${\rho }_m$ and consequent decrease of the stepped growth of ${\rho }_f$.

Now, consider the physical mechanism underlying the load drops. To understand this, we have calculated the plastic strain rate ${\dot{ϵ}}_p$ as a function depth. This is shown in Figure 5b. As can be seen, these bursts of $\dot{ϵ}$ occurs whenever ${\dot{ϵ}}_p$ exceeds $\frac{\dot{r}}{T}$, the imposed rate. It may also be noted that these bursts are well correlated with the bursts of ${\rho }_m$ as should be expected. The decreasing peak heights of ${\rho }_m$ or the load drops or ${\dot{ϵ}}_p$ is not due to the increasing back stress, but due to the depletion of ${\rho }_m$ as the indentation proceeds (see for details Ref. [28]).

Since our approach is based on dislocation mechanisms contributing to plastic flow, we can calculate the residual indentation depth by unloading the indenter from any point on the $F-z$ curve. When the indenter is unloaded beyond the first force drop, we find that the residual indentation depth, and consequently the residual imprint area is finite, another important feature of nanoindentation experiments. The red curve in Figure 4a shows the unloading curve of the indenter from the point marked • (after the fifth force drop). The residual imprint area is $\pi a^2$∼$18.78\times {10}^{-16}$ m${}^2$, a value comparable to the estimated residual contact area (from Figure 2a of Ref. [11]). Note however that the values of residual indentation depth $z_{pr}$ at the bottom of a load drop and the top of the following loading branch are nearly the same, but their load values are different. In contrast, $z_{pr}$ is different for points on the top and bottom of a load drop. Thus, assuming hardness is defined as the ratio of the load to the residual imprint area, hardness assumes two values and therefore is not well defined in the intermittent region. As a consequence, pile-up or sink-in effects are not relevant in the region of intermittent nanoindentation.

5.3. Summary and Discussion

In summary, we have demonstrated that the model not only predicts all the generic features of nanoindentation but also predicts numbers that match the experimentally measured load, stress, depth of indentation, the maximum stress on the Hertzian branch and residual plasticity under unloading. The good match between the numbers predicted by model and experiments can be attributed partly to the fact that our approach allows for a direct adoption of experimental parameters such as $\dot{r},R,T,b,E^{\ast }$ and $h=\alpha Gb$ and partly to the fact that we have included all the relevant dislocation mechanisms contributing to the intermittent indentation process. Indeed, the strength of our approach lies in providing an innovative method for calculating $z_p$ from the evolution equations for ${\rho }_m$ (and ${\rho }_f$) through the Orowan expression for the plastic strain rate.

Our model also provides insight into the physical and mathematical origin of the decreasing magnitudes of the load drops, a feature seen in experiments and in our model. This feature is not due to the back stress ${\sigma }_b=\alpha Gb{\rho }_f^{1/2}$ since the peak value of ${\rho }_f$ ranges from ∼5 to ∼7 $\times {10}^{13}$ m${}^{-2}$ giving ${\sigma }_b$∼40–47 MPa. At a physical level, this is purely due to the fact that every time there is a burst of ${\rho }_m$, much of it is transformed to ${\rho }_f$ depleting ${\rho }_m$ as indentation proceeds. Note also that there is no reverse transformation of ${\rho }_f$ to ${\rho }_m$. Such a reverse transformation is improbable due to the high stress required to break the junction. At a mathematical level, nanoindentation is a transient instability. The physical consequence of this is reflected in the decreasing magnitudes of bursts of ${\rho }_m$ and consequent load drops.

6. A Dislocation Dynamical Model for Displacement Jumps Using a Berkovich Indenter

For the LC mode indentation, we plan to address experimental results on Al single crystals using a Berkovich indenter. For this, we use appropriate expressions for the velocity of dislocation, the load, and the area specific to Berkovich indenter in Equations (15) and (16). For this case, we use the velocity expression given by Equation (14). Then, using this and Equations (10) and (11) in Equations (15) and (16), we obtain the evolution equations:

$$\begin{array}{ccc}\hfill {\dot{\rho }}_m& \hfill =& \frac{{\nu }_n}{\pi b^2}exp[\frac{z_{e_c}^2}{\frac{{\sigma }_n}{C^{\prime }{E}^{\ast }}[{\alpha }_1{(z_{e_c}+(z_p/s))}^2+{\alpha }_2(z_{e_c}+(z_p/s))]}]+\theta V_0{\rho }_m[\frac{\frac{C^{\prime }{E}^{\ast }{z}_{e_c}^2}{[{\alpha }_1{(z_{e_c}+z_p/s)}^2+{\alpha }_2(z_{e_c}+z_p/s)]}-h{\rho }_f^{1/2}}{{\sigma }_m}]{}^m\hfill \\ & \hfill -& {\nu }_{p_0}{\rho }_m[\frac{z_{e_c}^2}{\frac{{\sigma }_p}{C^{\prime }{E}^{\ast }}[{\alpha }_1{\left(z_{e_c}+z_p/s\right)}^2+{\alpha }_2(z_{e_c}+z_p/s)]}]{}^p-\beta {\rho }_m^2-\delta {\rho }_m{\rho }_f^{1/2}-f\beta {\rho }_m{\rho }_f,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\dot{\rho }}_f& \hfill =& \beta {\rho }_m^2-f\beta {\rho }_m{\rho }_f+\delta {\rho }_m{\rho }_f^{1/2}.\hfill \end{array}$$

(21)

Using the scale factor s between $z_e$ and $z_{e_c}$, Equation (12) takes the form

$$\begin{array}{ccc}\hfill \frac{d\left(C^{\prime }{E}^{\ast }{z_{e_c}}^2\right)}{dt}& =& {\dot{F}}_0.\hfill \end{array}$$

(22)

Note that we have assumed that the area-depth relation holds for the elasto-plastic region, i.e., $A\left(z\right)=A(z_{e_c}+z_{p_c})$, where $z_{p_c}$ is the contact depth contribution from plastic deformation.

The contribution arising from the plastic displacement $z_p$ is calculated by integrating

$$\begin{array}{ccc}\hfill {\dot{z}}_p& =& Tb{V}_0{\rho }_m[\frac{\frac{C^{\prime }{E}^{\ast }{z}_{e_c}^2}{[{\alpha }_1{(z_{e_c}+z_p/s)}^2+{\alpha }_2(z_{e_c}+z_p/s)]}-h{\rho }_f^{1/2}}{{\sigma }_m}]{}^m.\hfill \end{array}$$

(23)

Note that we have used $z_{p_c}=z_p/s$. Using the value of ${\nu }_s=0.33$ in Equation (10) and using the value of $s=1.6$ obtained by using Equation (8) of Ref. [38], we get $C^{\prime }=5.12704$.

Equations (20)–(23) constitute a coupled set of nonlinear ordinary differential equations for the LC mode nanoindentation problem. Such equations often exhibit instabilities for a domain of values of the model parameters. (See also Appendix of Ref. [28] where details of the stability analysis are given). In the present case, the instability domain of parameter values has been determined numerically. This shows that the instability occurs over a wide range of values of the parameters [29]. In general, the physically relevant values of the parameters form only a small subset of the instability domain. In the following, we briefly summarize the range of physically allowed values determined in Ref. [29].

6.1. Estimation of Parameters

As in the case of the DC mode indentation model [28], we adopt experimental parameters such as $E^{\ast },{\dot{F}}_0,R,b,T$ and $h=\alpha Gb$ and other shape parameters defining the Berkovich indenter geometry for single crystals of Aluminum [12]. The ranges of the parameters have been determined in our earlier publications [28,29] and is listed in

Table 2. Range of parameter values in the model equations using Berkovich indenter keeping the nominal tip radius $R=40$ nm and ${\sigma }_n=7.0$ GPa.

${\mathit{E}}^{\ast }$(GPa)	${\mathit{\sigma }}_{\mathit{m}}$(GPa)	${\mathit{\sigma }}_{\mathit{p}}$(GPa)	${\dot{\mathit{F}}}_0$($\mathsf{\mu }$N/s)
$73.7$	$0.7-1.0$	${\sigma }_m-6.5$	$0.3$
${\mathit{v}}_{{\mathit{p}}_{\mathbf{0}}}$ (m/s)	$\mathit{\beta }$ (m${}^{\mathbf{2}}$/s)	$\mathit{\delta }$ (m/s)	$\mathit{f}$
10${}^{-9}$–5 × 10${}^{-6}$	${10}^{-13}$	(0.01–4) × 10${}^{-6}$	${10}^{-2}$–10

. Very briefly, the nucleation stress ${\sigma }_n$∼$E^{\ast }/10$, the multiplication stress ${\sigma }_m$∼$E^{\ast }/300-E^{\ast }/100$, the propagation stress ${\sigma }_p<{\sigma }_m$. Parameter ranges of $\beta ,\delta$ and f have been determined along the lines mentioned for the DC mode in Section 5.1.

Recall that, for the LC mode, we have used the power law expression for the dislocation velocities given by Equation (14). This contains another parameter, namely, the velocity exponent m. This expression has been suggested on the basis of fits of dislocation velocities for bulk specimens subject to various loading conditions such as tensile and pulse loading using transmission electron microscopy, etch-pit technique, and slip line cinematography. (See Figure 19 of Ref. [48]). Under these conditions, stress is nearly uniform within the sample. Equation (14) is valid for single and groups of dislocations. Since the value of m depends on the material, the allowed range of values of m is difficult to estimate on physical grounds, particularly when conditions of deformation are very different (as for nanoindentation) from the conditions where it is measured. On the other hand, often, an inverse correspondence of the velocity exponent m with strain rate exponent n in $n=[\frac{dln\sigma }{dln\dot{ϵ}}]$ is suggested. For single crystals of Al, the value of n is ∼0.02 or less. This implies that the range of values of the velocity exponent m can be as high as 50. However, this correspondence breaks down in nanoindentation experiments at indentation depths ∼50 nm where displacement jumps dominate. Moreover, at such small length scales, it is well known that the yield stress increases rapidly, which is also suggestive of the fact that m values for bulk may not be valid for small indentation depth. In view of this, we have used m to range from 2–15.

Similarly, the exponent p (in the expression for velocity of propagation) is taken to be in the range from 1–10. However, p turns out to be an insensitive parameter. The range of values of the parameters are listed in Table 2.

Interestingly, it is possible to calculate the depth at which nucleation occurs given ${\sigma }_m$ and R the nominal radius of the blunted tip of the indenter or vice versa [29]. This is given by

$$z_{e_c}\left(t_n\right)\left[1-\frac{{\sigma }_n}{E^{\ast }}\frac{24.54}{s^{2}C}\right]=\frac{2\pi R}{s^{2}C}\frac{{\sigma }_n}{E^{\ast }}.$$

(24)

(Here, $s=1.6$.) Thus, the maximum elastic indentation depth $z_e\left(t_n\right)$ with nucleation occurring at time $t_n$ is completely determined by ${\sigma }_n/E^{\ast }$ and R or vice versa. Note that this statement is independent of the mode of deformation. Thus, Equation (24) allows determination of either ${\sigma }_n$ or $z_e\left(t_n\right)$ given the other assuming R is given by experiments. For example, if we use ${\sigma }_n/E^{\ast }=1/10$ and $R=40$ nm (taken from Ref. [12]), we get $z_n=14.9$ nm (see Table 2). Furthermore, given a value of ${\sigma }_n$, we can use $z_e\left(t_n\right)$ in Equation (10) to determine F. Then, Equation (12) determines the time $t=t_n$ at which the nucleation occurs. Thus, both F and $z_e\left(t_n\right)$ are also uniquely determined by Equation (24).

We can take the arguments further to show that the first displacement jump is controlled by ${\left[({\sigma }_n/{\sigma }_m)\right]}^m$. To do this, we note that dislocations are nucleated only when the stress $\sigma$ exceeds the nucleation stress ${\sigma }_n$. On the other hand, the magnitude of the displacement jump or load drop is directly controlled by ${[{\sigma }_n/{\sigma }_m]}^m{\rho }_m(t=t_n),$ where ${\rho }_m(t=t_n)$ is the dislocation density at the nucleation point. Note that we have made use of the fact at nucleation $\sigma ={\sigma }_n$ and that ${\sigma }_n$ is much larger than ${\sigma }_m$. Therefore, ${\left[({\sigma }_n/{\sigma }_m)\right]}^m$ determines the magnitude of the first displacement jump or vice versa. However, since ${\sigma }_n$ is already determined, ${\sigma }_m$ and m are completely determined by the magnitude of the first displacement jump. (See Ref. [29] for details.)

One interesting consequence of this is that the first displacement jump (or load drop) can only exist provided the elastic branch exists or nucleation should occur prior to the jump. This can be seen by noting that the factor ${[{\sigma }_n/{\sigma }_m]}^m{\rho }_m(t=t_n)$ that controls the magnitude of the first displacement jump decreases as we decrease ${\sigma }_n$ and eventually vanishes as ${\sigma }_n\to 0$. (Note that the elastic branch gets shorter and shorter is equivalent to ${\sigma }_n\to 0$.)

6.2. Results

Equations (20)–(23) are solved using adaptive step Runge–Kutta solver (MATLAB ‘ODE 15s’) with initial conditions ${\rho }_m=0$ and ${\rho }_f=0$ to mimic the absence of dislocations in the indented volume.

6.3. Influence of Parameter Variation on the Model Force–Displacement Curves

The two main objectives, namely, to build a model for the LC mode indentation that would predict all the generic features of the LC mode experiments and to examine the ability of the model to provide a good fit to a given experimental load–displacement data require that we first identify the dislocation mechanisms (or equivalently the corresponding rate constants) controlling the different regions of the load–displacement curve. Towards this end, we begin by investigating the influence of each of these parameters on the model $F-z$ curves. Recognizing that it is a multi-parameter space, it is convenient to fix the parameters at specific values (for instance, those listed in

Table 3. Optimized parameter values used for obtaining the best fit with the experimental load–displacement curve for single crystals of $\left(110\right)$ orientation [12]. Other parameters are fixed at ${\sigma }_n=7.0\mathrm{GPa},{\sigma }_m=0.91\mathrm{GPa}$ and $m=11$. Although the parameters ${\sigma }_p,v_{p_0}$ and p have little influence on the model $F-z$ curves, the values cited are the values used for numerical fit with the experiment.

${\mathit{\sigma }}_{\mathit{p}}$ (GPa)	${\mathit{v}}_{{\mathit{p}}_0}$ (m/s)	p	$\mathit{\beta }$ (m${}^2$/s)	$\mathit{\delta }$ (m/s)	f
$1.5$	$4\times {10}^{-7}$	4	${10}^{-13}$	$1.58\times {10}^{-6}$	10

) and study the influence of each of the parameters on the model $F-z$ curve. The study helps us to evaluate the relative importance of various dislocation mechanisms. In addition, the general trends of the influence of the parameters on the model $F-z$ curves would also be helpful to obtain an optimized set of parameters that provides a good fit to a given experimental $F-z$ curve.

Recall that the analysis reported in the previous section demonstrates that ${\sigma }_n$ can be taken to be fixed. Thus, we first consider the influence of ${\sigma }_m$ on the model $F-z$ curve with all other parameters fixed as in Table 3. (Note that ${\sigma }_m$ and m together specifies the dislocation multiplication mechanism. Therefore, while studying the influence of one, the other must be held fixed). Noting that, for small volume systems, ${\sigma }_m$ is larger than the bulk value, we have varied ${\sigma }_m$ from ∼$E^{\ast }/100$ to $E^{\ast }/70$. Plots of the $F-z$ curves for $E^{\ast }/{\sigma }_m=105,92,82$ and 74, are shown in Figure 6a. It is clear that increasing ${\sigma }_m$ decreases the magnitude of the first displacement jump and also the secondary jumps. More importantly, the model $F-z$ curves for all values of ${\sigma }_m$ exhibit displacement jumps of decreasing magnitudes, a characteristic feature of the LC mode indentation. Note that the load remains practically constant during the first displacement jump.

Now, consider the influence of the velocity exponent m on the model $F-z$ curve. Plots of $F-z$ curves for $m=4,8,10$ and 12 are shown in Figure 6b. The figure shows that the higher the value of m, the larger is the first displacement jump while the subsequent displacement jumps decrease marginally.

Note also that increasing ${\sigma }_m$ and increasing m have the opposite effect on the magnitude of the first displacement jump. This feature can be easily understood by noting that the first displacement jump is determined by ${[\frac{{\sigma }_n}{{\sigma }_m}]}^m$. Clearly, as ${\sigma }_m$ increases, the factor ${[\frac{{\sigma }_n}{{\sigma }_m}]}^m$ decreases while it increases with increasing m values. The opposing influence of ${\sigma }_m$ and m on the model $F-z$ curves turns out to be very helpful in obtaining the optimal values of ${\sigma }_m$ and m that match the magnitude of the first displacement jump for a given experimental data, in particular while fitting the experimental plots of single crystals of Al.

We have also examined the influence of $\delta$ (in the range $\delta =0.06$ to $3.16\times {10}^{-6}$ m/s ) for a range of f values in the interval ${10}^{-2}$ to 10. (Note that we can vary the product $f\beta$ or vary f keeping $\beta$ fixed). For small f, say 0.01 (or 0.1), as we vary $\delta$, the model $F-z$ curves exhibit several displacement jumps of decreasing magnitudes as can be seen in Figure 6c. Similarly, we have studied the influence of f on the model $F-z$ curve keeping $\delta$ fixed in the range $\delta =0.06$ to $3.16\times {10}^{-6}$ m/s . Plots of the $F-z$ curves for various value of f are shown in Figure 6d. Again, all model $F-z$ curves exhibit several displacement jumps of decreasing magnitudes. Furthermore, this study also shows that while larger values of $\delta$ lead to more rapid increase in the load after each displacement jump, larger f leads to larger displacement jumps with smaller load increase. Note the opposing trends of the influence of f and $\delta$ on the $F-z$ curves.

On the other hand, several other parameters ${\sigma }_p,p$ and ${\nu }_{p_0}$ have very little influence on the model $F-z$ curves. This is fortuitous considering the fact that these parameters were difficult to estimate.

In conclusion, the study of the influence of the parameters on the model $F-z$ curve demonstrates that for a wide range of values of the parameters, the model $F-z$ curves capture the characteristic features of the LC mode nanoindentation. In addition, the study determines the relative importance of different dislocation mechanisms, or equivalently, the relative importance of the corresponding parameters on the nanoindentation process.

6.4. Determination of Optimal Parameter Values for Best Fit to Experimental Load–Displacement Curve

The results of the previous section will now be used to obtain an optimized set of parameters that give best fit to the desired experimental load–displacement curve. The following four results from the previous section will be used for this purpose: (a) There is a range of values of the parameters in the instability domain for which the model $F-z$ curves exhibit all the characteristic features of experimental $F-z$ curves. This also implies that one should expect to find an optimized set of parameters within this range that fit the experimental data; (b) The study identifies three dislocation mechanisms contributing to the load–displacement curve in a sequential way. The maximum elastic depth $z_e\left(t_n\right)$ at which the nucleation occurs is determined by Equation (24). The nucleation stress ${\sigma }_n$ is determined by the maximum elastic depth $z_n$ or vice versa; (c) The magnitude of the first displacement jump is determined by ${[{\sigma }_n/{\sigma }_m]}^m$. However, since ${\sigma }_n$ is already determined, ${\sigma }_m$ and m are completely determined by the magnitude of the first displacement jump; and (d) The rest of the $F-z$ curve with multiple displacement jumps is determined by f and $\delta$ given the maximum experimental load and indentation depth.

We shall use this information to fit the load–displacement curve for single crystals of Al for $\left(110\right)$ orientation [12]. In experiments, the thickness of the samples used are 400, 600 and 1000 nm, and the radius of the indenters used are 50, 100, 150 nm. For the $\left(110\right)$ case, we have used $T=400$ nm and $R=40$ nm. To do this, we shall use the following data extracted from the experimental $F-z$ curve (Figure 7 of Ref. [12]) to determine the optimal parameters for the best fit. The maximum depth on the elastic branch is $z_e\left(t_n\right)\approx 14$ nm and the first displacement jump of $\Delta z_1\approx 33$ nm. The maximum load and maximum displacement are respectively $F_{max}\approx 44$ $\mathsf{\mu }$N and $z_{max}\approx 60$ nm. Using $z_e\left(t_n\right)\approx 14$ nm (and $R=40$ nm) in Equation (24) gives the nucleation stress ${\sigma }_n=7.0$ GPa. Now, consider finding the optimal values of ${\sigma }_m$ and m. Noting that increasing ${\sigma }_m$ decreases the first displacement jump and increasing m increases the magnitude of the first displacement jump, and using $\Delta z_1\sim 33$ nm gives ${\sigma }_m=0.91$ GPa and $m=11$. (Indeed, a careful perusal of the plots shown in Figure 6a,b also suggest value close to these values). These values give a good fit to the experimental $F-z$ curve until the end of the first displacement jump as is clear from Figure 7a. (Model $F-z$ curve is shown by the continuous line. Points marked • in the figure are experimental points extracted from Figure 7 of Ref. [12].)

We now consider getting a good fit to the rest of the experimental curve until $F_{max}\approx 44$ $\mathsf{\mu }$N and $z_{max}\approx 60$ nm. Recall that the portion of the $F-z$ curve beyond the first displacement jump is controlled by f and $\delta ,$ and therefore we need to find the optimized values of f and $\delta$ subject to the condition that the model $F_{max}$ and $z_{max}$ match closely the experimental values. To do this, we use the fact that $\delta$ and f have opposing influence on the model $F-z$ curve. Then, the values of f and $\delta$ that satisfy this condition are listed in Table 3.

Figure 7a shows the model $F-z$ curve (continuous curve) obtained by using the optimized parameter values along with the experimental points •. It is clear from the figure that the experimental points fall on the predicted curve except for the last few nanometers. The magnitudes of the successive displacement jumps decrease and so do the magnitudes of the load steps. Note also that the predicted load remains practically constant during the first displacement jump unlike in experiments where it decreases from 29.5 to 28.5 $\mathsf{\mu }$N [12]. The latter can be attributed to the inability to enforce constant load rate in experiments during the displacement jump.

Following the above procedure, we have also attempted to get the best fit to the experimental $F-z$ curve for Al single crystals of $\left(133\right)$ orientation. Both the model $F-z$ curve (continuous curve) and experimental points • are shown in Figure 7b. As can be seen, the fit is very good. Thus, the model captures not just the generic features of the experimental $F-z$ curve for a range of parameter values, the model $F-z$ curve corresponding to the optimized parameter values provide very good fit to the experimental data.

Our approach allows us to compute the stress as a function of time or indentation depth. This is shown in Figure 8a corresponding to the $\left(110\right)$ orientation. As can be seen, the maximum stress on the elastic branch $\sim 7.1$ GPa is close to the theoretical strength of the material. Thereafter, the stress relaxes largely during the first displacement jump. However, it must be noted that this relaxation occurs in a short duration of time, which also corresponds to the duration of the first displacement jump. Subsequently, the stress relaxation slows down eventually reaching the asymptotic value of ∼1 GPa. Note that this value is higher than the multiplication threshold stress ${\sigma }_m=0.91$ GPa.

As emphasized earlier, the first displacement jump beyond the elastic branch is triggered by the abrupt multiplication of the nucleated dislocations. Figure 8b shows the plot of the mobile density as a function of time. As can be seen, the first burst is large (${\rho }_m$∼$2.5\times {10}^{14}/{\mathrm{m}}^2$) and its duration is short. The magnitudes of successive bursts decrease rapidly to a value $2.3\times {10}^{11}{\mathrm{m}}^{-2}$ with their durations increasing as is clear from the inset of Figure 8b. These features are characteristic features of collective dislocation activity. The decreasing magnitudes of the ${\rho }_m$ bursts and increasing duration of these bursts clearly suggests that collective effects become weak and eventually disappear in the region where displacement jumps disappear. From a physical point of view, the decreasing magnitudes of the burst of ${\rho }_m$ is caused by the common loss term, $f\beta {\rho }_m{\rho }_f$ that depletes both ${\rho }_m$ and ${\rho }_f$ while the storage terms $\beta {\rho }_m^2$ and $\delta {\rho }_f^{1/2}{\rho }_m$ depletes ${\rho }_m$. The net effect is that both ${\rho }_m$ and ${\rho }_f$ decrease to their asymptotic values.

Since our approach is based on dislocation mechanisms, we can calculate the residual indentation depth by unloading the indenter. This also allows us to calculate the hardness as the ratio of the load to residual imprint area. However, experimental load-depth plots for Al are limited to depths where displacement jumps dominate and therefore hardness is not well defined. On the other hand, model calculation can be performed for indentation depths where displacement jumps disappear. Hardness in this regime is shown in Figure 9a. It is clear that hardness decreases as a function of depth. In experiments, small pile-up is visible in the micrograph shown in Figure 9 of Ref. [12]. However, the authors do not report hardness and therefore no comparison is possible. We emphasize that the decreasing trend of H is not altered for other values of s. Thus, our model predicts the indentation-size effect as well. The asymptotic values of $H=1.98$ GPa is almost twice the stress at large z (see Figure 8a).

Here, a few comments are in order on the existing models for hardness and our own radically different approach for calculating hardness. All traditional models of hardness are based on extending the Taylor relation for hardness to include the additional resistance arising from geometrically necessary dislocations (GNDs). The Nix–Gao model and its variants are based on Taylor relation, and therefore these models do not calculate hardness as a function of indentation depth, because, this amounts to calculating the density of geometrically necessary dislocations (GNDs) and statistically stored dislocations (SSDs) as a function of indentation depth. Indeed, these models use GNDs and SSDs (or equivalently asymptotic hardness and the slope of hardness as a function of indentation depth) as fitting parameters to experimental data.

We have recently extended our approach to calculate hardness as a function of indentation depth from small depths ∼nm to few $\mathsf{\mu }$m [51]. While the traditional models are based on Taylor relation for flow stress, a property that reflects the resistance to dislocation motion, our approach relies on calculating residual indentation depth, a property that is a measure of ease of dislocation motion. Our model includes GNDs in addition to mobile and forest dislocations (or equivalently SSDs). The model [51] predicts all the characteristic features of hardness such as the decreasing nature of hardness with increasing depths, square of the hardness scaling inversely with indentation depth larger than 150 nm and its deviation for smaller depths. In addition, we have obtained good fits to the two widely cited hardness data for strain hardened polycrystalline Cu [52] and single crystals of Silver [53]. To the best of our knowledge, this is the first indentation size effect model that calculates hardness as a function of indentation depth independently.

6.5. Physical and Mathematical Mechanisms for the Pop-in Events

One important feature of the LC mode nanoindentation is the existence of several displacement jumps of decreasing magnitudes. While the origin of displacement jumps is attributed to bursts of plasticity, explaining the decreasing magnitudes of the displacement jumps quantitatively has remained a puzzle. Clearly, this feature can be addressed within the scope of our model since the displacement jumps are controlled by bursts of ${\rho }_m$. Recall that studies on the influence of the parameters on the model $F-z$ curves (Section 6.3) has demonstrated that three different dislocation mechanisms operate in a time sequence. The initial elastic branch is controlled by nucleation of dislocations while the first displacement jump is controlled by the multiplication of the nucleated dislocations. On the other hand, the subsequent displacement jumps of decreasing magnitudes are controlled by the storage and recovery dislocation mechanisms. While the magnitudes of the displacement jumps are determined by the magnitudes of ${\rho }_m$ bursts, these themselves decrease. The latter is due to the fact that a large part of every ${\rho }_m$ burst is transformed to ${\rho }_f$. Furthermore, noting that the stress relaxes largely during the first displacement jump, the magnitudes of subsequent bursts of ${\rho }_m$ are determined by the excess stress above ${\sigma }_m$. Then, the decreasing over stress leads to decreasing magnitudes of the secondary ${\rho }_m$ bursts. This in turn implies decreasing magnitudes of the displacement jumps. This explanation should be contrasted with the one suggested in the literature, namely, increasing back stress. This is not the case since the back stress ${\sigma }_b=\alpha Gb{\rho }_f^{1/2}$ is ∼25 MPa even when ${\rho }_f$ assumes a peak value of ∼$3\times {10}^{13}/{\mathrm{m}}^2$.

6.6. Common Mechanism Underlying All Displacement Jumps

Recall that, in the DC mode of indentation, the load drops occur whenever the plastic displacement rate ${\dot{z}}_p$ exceeds the applied displacement rate $\dot{r}$ as is clear from the Equation (19). The equation explains all load drops on the same footing. However, an equivalent mechanism is absent for the case of the LC mode nanoindentation. In view of this, we attempt to explain all displacement jumps on the basis of a common dynamical mechanism. Some insight can be obtained by examining the stress plot in Figure 8a together with the $F-z$ plot in Figure 7a. Figure 8a shows that the stress relaxes mostly during the first pop-in event itself, eventually reaching its asymptotic value of ∼1 GPa. However, this value is larger than the dislocation multiplication stress ${\sigma }_m=0.91$ GPa. Since the stress remains higher than ${\sigma }_m$ even after the first displacement jump, one should expect a continuous increase in ${\rho }_m$ or equivalently a continuous increase in $z_p$. In contrast, Figure 7a shows that the load increases quite sharply beyond the first displacement jump suggesting a dominant elastic contribution. This raises a question as to why the load should increase in a near elastic manner when the stress remains higher than ${\sigma }_m$? As we shall see, this question is also relevant for subsequent displacement jumps.

To see this, recall that both load and area evolve with time. The area increase can occur either due to imposed loading rate or due to plastic displacement rate. However, the time scales of these two processes are very different. While the load increases linearly, the rate of area increase $\dot{A}$ arising from plastic deformation can be very short as is clear from the short time scale of ${\rho }_m$ bursts. From a dynamical point of view, this observation suggests that, while the underlying dislocation mechanisms for the first and secondary displacement jumps are different, we can describe all displacement jumps as resulting from a competition between the two well separated time scales, namely, the slow time scale corresponding to ${\dot{F}}_0$ and the fast time scale corresponding to ${\rho }_m$ bursts. The rate of area change arising from imposed load rate ${\dot{F}}_0$ is given $\frac{{\dot{F}}_0}{{\sigma }_m}$, where ${\sigma }_m$ is a natural choice for converting the applied load rate to the rate of area increase. Figure 9b shows a plot of $\dot{A}\left(z\right)$ as a function of time, which increases in bursts due to ${\rho }_m$ bursts. Also shown is the rate of area increase due to applied load rate $\frac{{\dot{F}}_0}{{\sigma }_m}=\mathrm{constant}$ (red line). It is clear that, whenever $\dot{A}\left(z\right)$ overshoots $\frac{{\dot{F}}_0}{{\sigma }_m}$, we see a displacement jump (triggered by bursts of ${\rho }_m$) that is almost fully plastic. Otherwise, the elastic component dominates over the plastic contribution. This feature is very similar to the mechanism underlying the force drops in the DC mode of indentation where a force drop occurs whenever the plastic displacement rate exceeds the applied displacement rate. Note that the duration of the first burst of $\dot{A}\left(z\right)$ (due to the first burst of ${\rho }_m$) is very short while the durations of the subsequent bursts get longer and longer until the magnitudes of $\dot{A}\left(z\right)$ and $\frac{{\dot{F}}_0}{{\sigma }_m}$ are comparable. Figure 9b can also be used to obtain quantitative estimates of the magnitudes of the load steps seen in Figure 7a. (See Ref. [29] for details.)

The above analysis demonstrates that the mathematical mechanism causing the instability is a result of a competition between the slow applied time scale corresponding to load rate (strictly $\frac{{\dot{F}}_0}{{\sigma }_m}$) and another fast time scale corresponding to $\dot{A}\left(z\right)$(controlled by bursts of ${\rho }_m$). The existence of one slow loading time scale and one fast time scale of response are typical of all relaxational oscillations (see Refs. [30,44,45]). However, in the present case, the instability is a transient one as is clear from the decreasing magnitudes of the displacement jumps and load steps. The analysis also demonstrates that imposing a constant load rate would be sufficient to trigger displacement jumps provided it is much larger than the short duration of ${\rho }_m$ bursts. However, in experiments, enforcing constant force rate condition during displacement jumps may not be easy due to finite response time of the machine.

6.7. Conclusions

In summary, the model for the LC mode nanoindentation predicts all the generic features such as the existence of an elastic branch followed by several displacement jumps of decreasing magnitudes, and residual plasticity under unloading for a range of values of the parameters. Furthermore, the predicted values of the load, the magnitudes of displacement jumps, etc. are also similar to those seen in experiments. The study of the influence of the parameters on the model $F-z$ curves show that the elastic depth is determined by ${\sigma }_n/E^{\ast }$ and the magnitude of the first displacement jump is controlled by the ratio of the nucleation stress to the multiplication stress, more precisely ${({\sigma }_m/{\sigma }_n)}^m$. Subsequent displacement jumps are controlled by the storage and recovery mechanisms. This identification also allows us to determine the optimized parameters that fit closely the experimental $F-z$ curves corresponding to single crystals of Al for $\left(110\right)$ and $\left(133\right)$ orientations.

Even though the underlying dislocation mechanisms for the first and subsequent jumps are different, we have demonstrated explicitly that all displacement jumps result whenever the short duration bursts of area increase due to plastic deformation exceed the slow rate of area increase due to applied load rate. Otherwise, elastic response dominates. Note that the short duration of bursts of ${\rho }_m$ is also a signature of collective behavior of dislocations. The mathematical mechanism for the instability is attributed to the competition between the slow time scale corresponding to the load rate and the fast time scale corresponding to bursts of mobile density. More specifically, the decreasing magnitudes of the displacement jumps or the load steps are signatures of transient nature of the instability. The underlying physical mechanism for decreasing magnitudes of the displacement jumps is the depletion of the mobile density that transforms to the forest density during the short duration of ${\rho }_m$ bursts.

7. General Conclusions on the Dynamical Approach to Nanoindentation

The present dynamical approach has been developed as an alternate approach to simulations with a view to overcome their inherent limitations. While these methods have been routinely used due to their ability to include a range of dislocation mechanisms [21,22,23,24,25,26], this advantage is offset by the serious limitations due to the short time scales inherent to several simulations and their inability to adopt experimental parameters, particularly the imposed rate of deformation. Consequently, the predicted numbers differ considerably from those reported by experiments. Our idea is to develop an alternate method to overcome the above limitations by devising a method that employs experimental parameters such as the rates, the geometrical parameters defining the geometry of the indenter, thickness of the film, etc., to predict numbers that match the experiments. The fact that the load drops or displacement jumps are signatures of an instability arising from collective dislocation activity motivated us to use nonlinear dynamics to model the nanoindentation instabilities. This also means that a natural interpretation of the dislocation densities used in the model is that they represent collective modes of the plastic deformation (as in the case of the AK model for the PLC effect).

The fact that the two models predict all the generic features with the predicted numbers matching those from experiments is clearly a support for both the dynamical approach and the dislocation mechanisms used. In this context, the success of the models in predicting all the generic features along with good fits to the two experimental data sets must be viewed against the background of the absence of simulations of any kind for the LC mode nanoindentation experiments. The fact that the model $F-z$ curves exhibit the characteristic features of the experiments for a range of parameter values can only be attributed to the fact that our approach allows for a direct adoption of the experimental parameters. More importantly, the strength of our approach lies in providing an innovative method for calculating the contribution from plastic deformation to the total indentation depth using the time evolution equations for ${\rho }_m$ and ${\rho }_f$. Note that obtaining a good fit to experiments requires that we match the model time and length scales with those in experiments. Noting that the total indentation depth is the sum of the elastic and plastic contributions, this match has been accomplished by demanding that the basic time scale for the growth of mobile density (that determines the extent of plastic deformation) matches the experimental time scale, i.e., we demand $\theta V_0$ is set to unity. (Note that imposed deformation such as the displacement rate or load rate or strain rate, are measured per second). This allows us to match the model $z_{max}$ with that in experiment. Note that matching model $F_{max}$ with that in experiment is straightforward since it is determined by load-elastic depth expression. Thus, our approach provides a consistent scheme to match the model time scale and length scale with those of the experiments. Note that the maximum depth is essentially determined by the plastic contribution to the depth.

One important consequence of our ability to calculate the plastic contribution to the indentation depth is that for larger depths where the load drops or displacement jumps disappear, the model predicts that hardness decreases with depth. Indeed, since our model uses dislocation density for the calculation of plastic contribution, our method provides a way to describe indentation size effect and constitutes an alternate approach to the strain gradient theories [54,55,56,57].

This kind of dislocation density based approach opens-up the possibility of studying nanoindentation on experimental length and time scales that cannot be achieved in finer scales simulation techniques. Furthermore, unlike simulations that require heavy computational resources, the model equations can be solved on a desktop computer.

The results predicted by the models for the DC mode and LC mode nanoindentation, in particular, the good fit to the two experimental load–displacement curves for the LC mode indentation where no simulations exist, clearly supports our view that the proposed method can be used as an alternate approach to simulation methods. This view is further supported by the success of the recent dislocation mechanism based dynamical model for calculating hardness as a function of indentation depth. The Indentation size effect model predicts all characteristic features of hardness in addition to providing good fits to the two widely used hardness data on cold worked polycrystalline Cu [52] and single crystals of Ag [53].

The approach may be criticized for ignoring the spatially inhomogeneous nature of the deformation. (Note that, although the indentation depth z is used as a dynamical variable, the model does not explicitly include spatial degrees of freedom.) Therefore, it might come as a surprise that, despite the absence of spatial degrees of freedom in the model, the model $F-z$ curve matches the experimental data very well. To appreciate this, it is important to recognize that dislocation densities used here represent the spatial/sample averages over the indented volume. They also represent collective dislocation modes that are responsible for load drops or displacement jumps. Much the same way, all experimentally measured quantities such as the force, stress, depth of indentation, residual plasticity under unloading, etc. are also the volume averages of the dislocation activity in the sample. Since both theoretically computed and experimentally measured quantities represent sample averages, the good match is not all that surprising. Moreover, it is well known that averages are quite insensitive to the details of the distribution. This statement is applicable to spatial averages as well. However, if spatial degrees of freedom are introduced, displacement jumps will exhibit some stochasticity, i.e., displacement jumps as well as the steps on force will be different in different runs. This is a general result in dynamical systems that has been well illustrated in the case of the PLC effect (see Ref. [30,33] and in particular Figure 2.1a of Ref. [42]). Inclusion of spatial degrees of freedom for the present problem is nontrivial since this involves a moving boundary that goes by the name Stefan’s problem in mathematics. Finally, it would be interesting to explore the possibility of using the current approach along with finite element methods that calculate the stress distribution under the indenter. This then would allow the possibility of obtaining spatial distribution of the dislocation activity in the sample.