# Revisiting the statistical estimation of activation parameters in incipient
plasticity

Sweta Kumari  Amlan Dutta<EMAIL_ADDRESS>Department of
Metallurgical and Materials Engineering, Indian Institute of Technology
Kharagpur, West Bengal 721302, India

###### Abstract

Incipient plasticity is typically associated with thermally activated events
like the nucleation of dislocations in crystalline solids and the activation
of shear transformation zones in metallic glasses. A method of estimating the
activation parameters of such mechanisms is to analyze the statistical
distribution of critical loads obtained through a series of repeated
measurements. However, this approach has been observed to produce activation
volumes of the order of atomic volumes in a variety of materials and
experimental setups. Such exceptionally small activation volumes have been
explained by conjecturing a non-trivial mechanism of nucleation. Here, we
critically analyze the inherent assumptions of the statistical method and show
that unexpected activation volumes can emerge simply from the statistical
fluctuations in the activation parameters themselves. To this end, we perform
repeated deformation simulations of iron nanopillars under both tensile and
compressive loading and measure the resulting yield stresses. Although the
conventional statistical analysis exhibits extremely small atomic volumes, the
atomistic simulations indicate a transition pathway that is physically
incommensurate with the statistical result. Using a simple Monte Carlo scheme
and analytical consideration, we show that even a relatively small dispersion
in activation parameters can misleadingly suppress the measured activation
volume to a significant extent. This shows that the ultra-small atomic volumes
reported in the earlier studies do not need exotic mechanisms but can be
explained simply as the misleading result obtained by ignoring the physically
plausible reality of statistical dispersion of activation parameters.

##  I Introduction

Experiments like tensile or compressive loading of solids and nano-indentation
exhibit abrupt transitions from elastic to plastic modes of deformation in the
form of yield points or pop-in events [1, 2, 3, 4]. Such transitions are
typically facilitated by load-driven activated processes like nucleation of
dislocation loops [5, 6, 7], depinning of line defects [8, 9], formation of
shear-transformation zones in metallic glasses [10, 11], etc. Over a broad
variety of solids, the onset of plasticity is determined by a combination of
temperature and rate of deformation. This indicates that the underlying
process involves load-driven, thermally activated events. The role of
temperature is to provide the thermal energy required by the system to proceed
along the minimum energy path, whereas the applied load effectively reduces
the activation barrier along the transition pathway.

In the context of any generic load-driven process, the activation volume is a
fundamental parameter associated with the activated event. It represents the
rate at which the activation barrier reduces with an increase in the applied
load. However, this key parameter plays a broader role by offering insights
into the possible mechanism underlying an activated process, as the activation
volume also provides an intuitive idea about the number of atoms associated
with the event. While the direct measurement of activation volume is
straightforward in atomistic computations with methods like nudged elastic
band (NEB) [12, 13, 14], umbrella sampling [15, 16], finite temperature string
[17], etc., experimental methods resort to mathematical models linking
activation volume to other measurable parameters. For instance, the activation
volume of solid state diffusion has been estimated by measuring diffusivities
under varying pressure [18]. In metallic glasses, the activation volume for
STZ activities is indirectly determined by measuring the hardness and strain-
rate sensitivity [19]. Similarly, Kato [20] has demonstrated a line tension
model that can be used to estimate the activation volume of dislocation
depinning at grain boundaries using experimentally gathered flow stress data
at various temperatures.

In 2004, Schuh and Lund developed an interesting model to extract the
activation parameters of a load-driven process by analyzing the statistical
distribution of the critical loads of the elastic-plastic transition, obtained
through a series of identical experimental measurements [21]. While the
original method was derived and demonstrated for nano-indentation of
crystalline samples, the fundamental idea is quite generic and easily
applicable to a variety of materials and processes. As a consequence, this
statistical method has witnessed numerous applications in various studies
involving experiments and simulations. However, despite the widespread use of
this technique, the activation volumes rendered by this method have been
somewhat perplexing. A number of studies, mostly using nano-indentation
experiments, have estimated activation volumes on the order of atomic volumes.
Such activation volumes are too small to be expected in dislocation mediated
processes. Several possible scenarios have been proposed to justify these
strange results. In the original study, Schuh and Lund suggested the
possibility of having nuclei of sub-equilibrium size. Subsequently, Mason et
al. expressed the possibility of the presence of sub-critical dislocation
loops and the effect of image stresses. Irrespective of the variety of
explanations, the statistical method has consistently produced extremely small
activation volumes in both crystalline and amorphous solids.

One common trend in all the studies based upon the statistical method of
obtaining the activation parameters has been that the observation of
unexpectedly small activation volumes has either been ignored or explained by
introducing the possibility of some non-trivial effect involved in the process
of deformation. To the best of our knowledge, none of the existing studies has
critically analyzed the effect of the various assumptions underlying the
statistical model. Mason et al.[22] considered the possible effects of surface
defects that could affect the resulting statistics of pop-in loads measured
during the nano-indentations. However, such factors are implicitly assumed to
have only secondary effects on the statistical distribution, and the measured
parameters are expected to be estimated with reasonable accuracy.

In the present study, we revisit the fundamental premise of the statistical
method of determining activation energy and volume. In particular, we examine
the assumption that during the measurement of critical loads, the exact same
activated event with the same activation barrier occurs every time the
critical load is recorded. To this end, we perform atomistic simulations of
the tensile and compressive deformation of iron nanopillars. This system
offers an interesting case study on account of the tension-compression
asymmetry found in bcc crystals. Molecular dynamics (MD) simulations have
established that while the tensile deformation of the Fe nanopillar is
governed by twinning, the compressive deformation involves the nucleation and
glide of screw dislocations. A series of computational experiments produce the
datasets of yield loads, which are statistically analyzed, and the calculated
activation volumes are directly compared to the mechanisms observed during the
simulations.

##  II Atomistic Simulations

Large-scale MD simulations of uniaxial tension and compression tests are
performed on virtual Fe nanopillars. The inter-atomic forces between Fe atoms
are described by the embedded atom method (EAM) potential developed by
Mendelev and co-workers [23]. In addition to replicating the basic material
properties, this potential correctly reproduces the non-degenerate core-
structure of screw dislocations and displaced twin boundaries in agreement
with the DFT calculations.

Figure 1: The bcc-Fe nanopillar sample

The 15.4 nm long nanopillar has a 5.8 nm $\times$ 5.8 nm square cross-section
in between two larger indenter plates, as seen in Fig. 1. The <100>-oriented
pillar has low energy {110}-surfaces, which prevent some spurious deformation
effects associated with the high energy surfaces. Uniaxial deformation is
applied at a rate of $10^{7}\,\text{sec}^{-1}$ and 300K temperature with a
time-step of 2 fs for the velocity-Verlet integrator. A set of 80 simulations
with random velocity initializations have been carried out, each with tensile
and compressive loading. During the deformation, the forces on the indenter
plates are recorded with yield the stresses on the nanopillar.

In addition to the MD simulations, the mechanism of yielding has been examined
by means of nudged-elastic band calculations. To this end, the structures of
the system just before and after the yield point are extracted. These
structures are axially strained until the desired stress is obtained and are
then used to initialize the chain of replicas. This chain is allowed to evolve
using the modified inertial relaxation method with the climbing replica
scheme. The NEB computations reveal the pathways of the elastic-to-plastic
transitions in both tensile and compressive deformations. All the computations
reported in this work have been carried out using the LAMMPS MD code [24]. The
OVITO program [25] with dislocation extraction analyzer aids the visualization
and analysis of crystalline defects.

##  III Statistical Model

The statistical model developed by Schuh and Lund offers the advantages of
simplicity and an exact solution. The model is based on the realization that
for a given loading rate, the rate of increase of the cumulative probability
(_F_) of a transition event is given by the product of the survival
probability (1-_F_) and transition rate ($\nu$). Assuming that the activation-
free energy reduces linearly with increased load over the range of interest,
we have the transition rate,

$\nu=\nu_{0}e^{-\Delta G/kT},$ (1)

with $\nu_{0}$ as the attempt frequency. The free energy barrier, $\Delta
G=G_{0}-v^{\ast}\sigma$ depends on the activation volume, $v^{\ast}$, and the
relevant stress component, $\sigma$. Solving the resulting differential
equation for the cumulative probability renders the solution,

$F(\sigma)=1-exp\left[-\frac{\eta
kT}{\dot{\sigma}v^{\ast}}\exp\left(\frac{\sigma v^{\ast}}{kT}\right)\right],$
(2)

which can rewritten as,

$\ln{\left[\ln(1-F)^{-1}\right]\ =\ \ln{\frac{\eta
kT}{\dot{\sigma}v^{\ast}}+\frac{\sigma v^{\ast}}{kT}}},$ (3)

where $\eta=\nu_{0}e^{-G_{0}/kT}$ and $\dot{\sigma}=d\sigma/dt$. The basic
essence of the statistical model can be captured in a simple scheme of Monte
Carlo calculations. We assume that the load $\sigma$ increases iteratively in
a small step of 0.01 per iteration, starting with $\sigma=0$. Each iteration
corresponds to a time span of unity. In each iteration, the barrier $\Delta G$
is evaluated and a uniformly distributed random number is generated. If this
number is found to be less than $e^{-\Delta G/kT}$, the iterations are
terminated. In the context of Eq. (3), this scheme simulates a transition
event with $\nu_{0}$ = 1 and $\dot{\sigma}=0.01$ in arbitrary units. We
perform calculations with the parameters, $G_{0}\ =\ 1$, and _kT_ = 0.05. For
each selected value of the activation volume, the Monte Carlo calculations are
executed 500 times with different initializations of the random number
generator.

Fig. 2 displays the variation of the reciprocal of the complementary
cumulative distribution scaled to a double logarithmic form. For three
different activation volumes, $v^{\ast}$ = 0.8, 1.0, and 1.2, the plots
exhibit linear behavior as predicted by Eq. (3). The excellent agreement
between the distributions given by the Monte Carlo calculations and those
obtained from Eq. (3) verifies the mathematical correctness of the statistical
model and the validity of the Monte Carlo scheme.

Figure 2: The reciprocal of the complementary cumulative distribution is
scaled to a double logarithmic form of shear stresses that are associated with
the critical load.

##  IV Results and Discussion

###  IV.1 Deformation simulations of Fe nanopillars

Fig. 3 exhibits the typical stress vs. strain plots of the nanopillars under
tensile and compressive loadings. The transition from elastic-to-plastic
deformation is easily identified as the sudden drop in measured load following
the yield point. As observed in earlier studies, the plasticity of tensile
deformation is governed by the formation and growth of a twinned region. In
contrast, the compressive yielding is controlled by the nucleation and glide
of screw dislocations. The snapshots of simulation cells accompanying the
stress-strain plots in Fig. 3 reveal these two distinct mechanisms of
deformation.

Figure 3: The stress on the Fe nanopillar is plotted against both tensile and
compressive strain. The black plot represents the tensile loading, while the
red plot represents the deformation under compression.

As specified in Sec. II, multiple randomized simulations have been performed
for both types of deformation. Fig. 4 shows the cumulative probabilities (_F_)
of the yield stresses plotted as double-logarithmic scale similar to those in
Fig. 2.

Figure 4: The cumulative probabilities (_F_) of the yield stresses are plotted
on a double-logarithmic scale for both tensile and compressive stress. An open
circle represents tensile loading, while a closed circle represents
deformation under compression.

We note that the linear fit is valid only for a part of the plot with tails
deviating from the relation predicted in Eq. (3), though the results of
compressive loading provide better fits over a broader range than their
tensile counterparts. This tail effect observed in these MD simulations is
consistent with the results observed in nano-indentation experiments.The tails
are either described as _less reliable_ part of the measured cumulative
distribution or often omitted while presenting the results. Table I presents
the activation volumes in the units of atomic volume $(\Omega)$. It must be
pointed out that in Eq. (3), we directly use the normal stress, $\sigma$,
instead of its resolved shear component, $\tau$. As a result, the activation
volume referred to in this work is assumed to be defined as,
$v^{*}=-\frac{\partial G}{\partial\sigma}$, which is typically expected to be
2-3 times smaller than the conventional activation volume, $-\frac{\partial
G}{\partial\tau}$, on account of the Schmid-factor. We observe that the
activation volumes obtained from the compressive deformation are substantially
smaller than those obtained from the tensile deformation. This is apparently
in consonance with our intuitive understanding, as we expect the heterogeneous
nucleation of a full dislocation loop to involve fewer atoms as compared to a
twin nucleus involving multiple slip planes and a larger volume. However,
despite the fact that the statistical analysis succeeds in revealing the
tension-compression asymmetry, the perplexingly small activation volumes still
remain unexplained.

Table 1: The values of fitting parameters for tensile and compressive loading Mode of Deformation | $\eta(\text{m}^{-3}\text{s}^{-1})$ | $v^{*}\,(\Omega)$
---|---|---
Tension | $4.93\times 10^{5}$ | 3.09
Compression | $2.276\times 10^{16}$ | 1.834

Nudged-elastic band calculations provide a clear view of the activation
process involved in yielding under tensile and compressive loads. As observed
in Fig. 5(a), the nucleation of twin fault takes place at the surface, and the
twinned region grows gradually in tensile deformation. Similarly, in
compressive deformation Fig. 5(b), the dislocation loop is seen to nucleate at
the corner, which soon re-orients to a screw-character due to the drastic
difference in mobilities of edge and screw components [26]. The transition
pathways for both types of deformation indicate that the processes are
mediated by planar slips, as per expectations. Furthermore, the saddle states
involve large enough volumes of the material involved in the activation
process, which is antithetical to the unusually small activation volume
observed in Table I.

The original statistical framework for estimating the activation parameters
was primarily employed for experimental observations. As it has not been
possible to actually witness the underlying mechanism of plastic yielding,
some non-trivial mechanisms have been proposed to explain the ultra-small
activation volumes. One of the possibilities that were expressed pertains to
the nucleation of numerous point defects beneath a nano-indenter tip. Another
possibility raised by the authors was the nucleation of extremely small loops
with activation volumes of the order of atomic volumes. In the subsequent
study, the mechanisms involving the migration of point defects were deemed
unfeasible. Instead, the surface asperities were proposed to activate the
heterogeneous nucleation of dislocation loops. Heterogeneous nucleation in a
supersaturated vacancy cloud was also discussed as a plausible route for pop-
in events. Zhu et al.[27] estimated the activation volume with first pop-in
statistics in fcc high-entropy alloys and found it to be around 3 atomic
volumes. They proposed that the incipient plasticity is governed by the
heterogeneous nucleation of dislocation loops with a critical size of only
$\sim 3$ atoms. A similar activation volume was reported for the bcc high-
entropy alloy [28] and the authors suggested the mechanism of cooperative
migration of a few atoms to be at play.

In spite of the exotic yielding mechanisms put forward to justify the
statistically estimated activation volumes, a detailed survey of the reported
results presents a perplexing observation. The statistical model has been
applied to a variety of experimental results obtained for pure metals,
ceramics, high-entropy alloys, and even metallic glasses. In each of these
studies, the activation volumes are invariably found to be of atomic order,
irrespective of the choice of material and mechanism of deformation. The same
trend is maintained in the reported MD simulations of Mo nanoparticles [29]
and of the Fe nanopillar as seen in the present work. As a matter of fact,
none of the reported studies has ever estimated a statistically computed
activation volume that is not of the order of atomic volume. Interestingly,
such activation volumes are always consistent with other methods of
measurement. For instance, Brandl et al. [30] studied the temperature-
controlled nano-indentation of bcc chromium and, using a different technique,
estimated the activation volume to be one or two orders of magnitude larger
than those typically obtained from the statistical method. Surprisingly, this
estimation was supposed to be for kink nucleation, which is supposed to
involve only a few atoms. Even if one ignores the loose interpretation of
activation volume in terms of the number of atoms involved in a critical
event, the straightforward definition still indicates that a very small
activation volume implies extreme insensitivity of the activation barrier
towards changes in the applied load. It is indeed confusing, as a dislocation-
mediated mechanism is supposed to be load-dependent and presents a strong case
in favor of reexamining the fundamental premises of the statistical approach.

###  IV.2 Revisiting the statistical model

The observations highlighted at the end of the previous subsection, along with
the activation paths revealed in Fig. 5, suggest that the exceptionally small
activation volumes estimated from the statistical analyses could rather be a
universal feature than the outcome of unusual mechanisms. This leads us to
reconsider the fundamental structure of the original statistical approach. As
the mathematical result is rigorous within the validity of its assumptions and
already verified by the Monte Carlo calculation (Sec. III), it is pertinent to
assess the fundamental premise of this model. We note that the model
intrinsically attributes the statistical variance of the measured critical
load to the thermally controlled probability of success of a barrier-crossing
transition, while the activation parameters remain the same for each
measurement. This can only be justified under one of the two possible
premises: first, if we assume that the exact same deformation process occurs
repeatedly during all the measurements; and second, if multiple processes are
deemed possible, albeit with some variation in activation parameters, the
statistical analysis is still assumed to yield a reasonable estimate of their
averages or some other representative values.

The validity of the first premise is apparently unlikely in a realistic
scenario. In a region of interest, _e.g.,_ beneath the tip of a nano-indenter,
there can be multiple sites where dislocation nucleation may occur with
different activation volumes owing to physical or compositional heterogeneity.
Similarly, nucleation on different atomic planes can occur with different
activation parameters. Another possibility is simultaneous nucleation at more
than one site, followed by interactions between the nucleated loops. Recent
studies have indicated that including more than one event type can affect the
statistical distribution. Tan and coworkers [31] have shown that the
generalized stacking fault energy in a complex concentrated alloy may be
modeled as smoothly transitioning between two values depending on the size of
the nucleated dislocation loop. This approximation can capture the effect of
spatial compositional fluctuation and provide a better fit to cumulative
distribution. Similarly, Nisanay and Mordehai [32] have theoretically analyzed
the situation of having multiple sites of dislocation nucleation. Their MD
simulations examined the nucleation of dislocation loops in nanowires with a
rhombic cross-section, thereby providing acute and obtuse sites for possible
heterogeneous nucleation. Even if we consider a single event type in a series
of experiments, it is possible to have more than one transition pathway
joining the initial and final states, with each path having its own saddle
point, activation barrier, and activation volume. In such a situation, the
total rate of transition would be the sum of the individual rates along the
transition pathways.

Figure 5: Nudged Elastic band images for (a)tension and (b)compression

As the first premise of having rigidly fixed activation parameters with no
statistical dispersion has been argued to be physically unrealistic, we can
consider the second possible premise. Here, we presume that the statistical
sampling of critical loads involves some inherent distribution of the
activation parameters. If the dispersion in the parameters is not too large,
the statistical approach can be considered to produce the correct parameter
within a sufficiently narrow range of tolerance. _Prima facie_ , this
assumption seems reasonable. Even if the repeated measurement of critical load
occurs under the possibility of multiple varieties of events, such events are
expected to be of the same fundamental nature. For instance, as observed in
our simulations, the tensile yielding of the nanopillar is always governed by
deformation twinning, and the compressive yielding proceeds via the nucleation
of isolated screw-type line defects, although the exact transition pathways
may have little variation. Similarly, the observation of first pop-in event
for a given material under nano-indentation is expected to occur through a
similar mechanism during multiple measurements. Thus, the activation
parameters can be considered not to vary in a drastic manner from one
measurement to the next.

The above discussion leads us to verify whether some dispersion in the
activation parameters within their given ranges can indeed produce an
acceptable statistical estimate of their representative values. To this end,
we may revisit the original statistical approach and introduce the statistical
variance of activation parameters into the model. Let us assume that a set of
experiments involves measurements of the critical load at a given temperature
and loading rate. Every measurement involves a critical event, _i_ , which is
picked from a pool of multiple possible events with probability,
$p(G_{0}^{i},\nu_{0}^{i},v_{i}^{\ast},\ \sigma)$, with the corresponding
activation parameters, $G_{0}^{i}$, $\nu_{0}^{i}$, and $v_{i}^{\ast}$(Sec.
III). If $\delta N_{i}$ is the number of measured load samples in the
infinitesimal range of $\sigma$ to $\sigma+d\sigma$, corresponding to event
_i_ , we have,

$dN_{i}=N_{i}f_{i}(\sigma)d\sigma,$ (4)

where $f_{i}(\sigma)$ is the probability density of the event of type, _i_ ,
with respect to variation in the measured critical load, $\sigma$. Summing
over all the possible event types, we obtain,

$dN=\sum_{i}{dN}_{i}=\sum_{i}{N_{i}f_{i}(\sigma)d\sigma}.$ (5)

Since the experimental or _in_ _silico_ load measurement cannot identify the
event type for a particular measurement, any direct estimate of the
probability distribution, $f_{i}$, is inaccessible. However, the measurements
still provide us with the distribution, _f_($\sigma$), of all the sampled
loads collectively and $dN=Nf(\sigma)d\sigma$. As,
$N_{i}=Np(G_{0}^{i},\nu_{0}^{i},v_{i}^{\ast},\ \sigma)$, Eq. (5) yields,

$f(\sigma)d\sigma=\sum_{i}{p(G_{0}^{i},\nu_{0}^{i},v_{i}^{\ast},\
\sigma)f_{i}d\sigma},$ (6)

and the cumulative distribution is obtained as,

$\Phi(\sigma)=\sum_{i}{\int_{-\infty}^{\sigma}{p(G_{0}^{i},\nu_{0}^{i},v_{i}^{\ast},\
\sigma)f_{i}(\sigma)d\sigma}}.$ (7)

It must be noted that the distribution,$\ \Phi(\sigma)$, includes all the
possible event-types, _i_ , and is fundamentally distinct from cumulative
distribution, $F(\sigma)$, which inherently assumes a single set of fixed
activation parameters. Although Eq. (7) is derived based upon the assumption
of a discrete set of possible parameters, {$\eta_{i},v_{i}^{\ast}$}, one can
easily generalize it to include a continuous distribution over the parametric
space as given below.

$\displaystyle\Phi(\sigma)$ $\displaystyle=$
$\displaystyle\int_{G_{0}^{\min}}^{G_{0}^{\max}}\int_{\nu_{0}^{\min}}^{\nu_{0}^{\max}}\int_{v_{\min}^{\ast}}^{v_{\max}^{\ast}}\int_{-\infty}^{\sigma}p(\eta,v^{\ast},\sigma)$
(8) $\displaystyle
f(G_{0},\nu_{0},v^{\ast},\sigma)\,d\sigma\,dv^{\ast}\,d\nu_{0}\,dG_{0}$

Equation (8) represents an extension of the original statistical scheme to its
most generic form. Various assumptions can be made to further reduce it to
simplified forms. For instance, if the probability distribution, _p,_ is
considered to be independent of the measured load, $\sigma$, Eq. (8) becomes,

$\displaystyle\Phi(\sigma)$ $\displaystyle=$
$\displaystyle\int_{G_{0}^{\min}}^{G_{0}^{\max}}\int_{\nu_{0}^{\min}}^{\nu_{0}^{\max}}\int_{v_{\min}^{\ast}}^{v_{\max}^{\ast}}p(G_{0},\nu_{0},v^{\ast})$
(9)
$\displaystyle\left\\{\int_{-\infty}^{\sigma}f(G_{0},\nu_{0},v^{\ast},\sigma)d\sigma\right\\}dv^{\ast}d\nu_{0}dG_{0}$

Here it must be pointed out that while Eq. (9) provides the cumulative
distribution over the whole 3-dimensional parametric space, Eq. (2) still
represents the partial cumulative distribution for a specific event-type with
fixed activation parameters, $G_{0}$, $\nu_{0}$, and $v^{\ast}$. Hence, Eq.
(9) can be rewritten as,

$\displaystyle\Phi(\sigma)$ $\displaystyle=$
$\displaystyle\int_{G_{0}^{\min}}^{G_{0}^{\max}}\int_{\nu_{0}^{\min}}^{\nu_{0}^{\max}}\int_{v_{\min}^{\ast}}^{v_{\max}^{\ast}}p(G_{0},\nu_{0},v^{\ast})$
(10)
$\displaystyle\left[1-\exp\left\\{-\frac{\eta_{0}kT}{\dot{\sigma}v^{\ast}}\exp\left(\frac{\sigma
v^{\ast}}{kT}\right)\right\\}\right]$ $\displaystyle dv^{\ast}d\nu_{0}dG_{0}.$

It is pertinent to note that the special case of a single event-type with
constant activation parameters, say $G_{0}^{c}$, $\nu_{0}^{c}$, and
$v_{c}^{\ast}$, is mathematically reproduced by setting the probability
density over the parametric space through Dirac-delta functions as,
$p(G_{0},\nu_{0},v^{\ast})=\delta(G_{0}-G_{0}^{c})\delta(\nu_{0}\nu_{0}^{c})\delta(v^{\ast}-v_{c}^{\ast})$,
thereby retrieving the original result of Eq. (2) as,

$\Phi(\sigma)=1-exp\left\\{-\frac{\eta_{c}kT}{\dot{\sigma}v_{c}^{\ast}}\exp\left(\frac{\sigma
v_{c}^{\ast}}{kT}\right)\right\\},$ (11)

with $\eta_{c}=\nu_{0}^{c}e^{-G_{0}/kT}$.

###  IV.3 Validation of extended model

The cumulative distribution expressed through Eq. (8) or even its moderately
simplified version in Eq. (10) still fails to enable the direct calculation in
the absence of _a priory_ knowledge of the distribution,
$p(G_{0},\nu_{0},v^{\ast})$, of the activation parameters. However, its effect
on the overall cumulative distribution can be observed by employing a simple
form, which can be easily verified through the Monte-Carlo scheme presented
earlier in Sec. III. The simplest way to introduce a variation in the
activation parameter is to consider the parameters to be uniformly distributed
around their corresponding mean values. We further assume that all the event
types have the same attempt frequency, which means that the probability
distribution of parameters only depends on the barrier-intercept, $G_{0}$, and
the activation volume, $v^{\ast}$. Thus, $\eta$ and $v^{\ast}$ are randomly
picked from the ranges $[\left\langle
G_{0}\right\rangle-\mathrm{\Delta}G_{0},\ \left\langle
G_{0}\right\rangle+\mathrm{\Delta}G_{0}]$ and ${[\left\langle
v^{\ast}\right\rangle-\mathrm{\Delta}v}^{\ast},\left\langle
v^{\ast}\right\rangle+\mathrm{\Delta}v^{\ast}]$, respectively, and
$p(G_{0},v^{\ast})=\frac{1}{4\mathrm{\Delta}G_{0}{\mathrm{\Delta}v}^{\ast}}$.
Then, the cumulative probability is extracted from Eq. (10) as,

$\displaystyle\Phi(\sigma)$ $\displaystyle=$ $\displaystyle\frac{1}{4\Delta
G_{0}\Delta v^{\ast}}\int_{\langle v^{\ast}\rangle-\Delta v^{\ast}}^{\langle
v^{\ast}\rangle+\Delta v^{\ast}}\int_{\langle G_{0}\rangle-\Delta
G_{0}}^{\langle G_{0}\rangle+\Delta G_{0}}$ (12)
$\displaystyle\left[1-\exp\left\\{-\frac{\eta
kT}{\dot{\sigma}v^{\ast}}\exp\left(\frac{\sigma
v^{\ast}}{kT}\right)\right\\}\right]\,$ $\displaystyle dG_{0}\,dv^{\ast}$

Equation (12) can be validated simply through small modification of the Monte-
Carlo method. At the beginning of each iteration, instead of using fixed
values of $G_{0}$ and $v^{\ast}$, we pick them randomly from a narrow range
around the central values, $\left\langle G_{0}\right\rangle$ and $\left\langle
v^{\ast}\right\rangle$. Fig. 6(a),(b) and (c) displays the of cumulative
distribution for $\left\langle G_{0}\right\rangle$ = 1 with
$\mathrm{\Delta}G_{0}$ as 5% of $\left\langle G_{0}\right\rangle$. Three
different central values of the activation volume, $\left\langle
v^{\ast}\right\rangle$= 0.8, 1.0, and 1.2 have been used with variations of
±5%, ±10%, and ±20%. Similarly, the results given in Fig. 6(d), (e) and (f)
correspond to ±5%, ±10%, and ±20% variation around $\left\langle
G_{0}\right\rangle$ with central values of $\left\langle
G_{0}\right\rangle$=0.8, 1.0, and 1.2 while the variation around $\left\langle
v^{\ast}\right\rangle$ is fixed at ±5% of $\left\langle
v^{\ast}\right\rangle$=1.

Figure 6: The of cumulative distribution for $\left\langle G_{0}\right\rangle$
= 1 with $\mathrm{\Delta}G_{0}$ as 5% of $\left\langle G_{0}\right\rangle$.
Three different central values of the activation volume, $\left\langle
v^{\ast}\right\rangle$ is used with variations of ±5%, ±10%, and ±20% (a)
$\left\langle v^{\ast}\right\rangle$= 0.8 (b) $\left\langle
v^{\ast}\right\rangle$= 1.0, and (c) $\left\langle v^{\ast}\right\rangle$=
1.2. The results (d), (e) and (f) correspond to ±5%, ±10%, and ±20% variation
around $\left\langle G_{0}\right\rangle$ with central values of $\left\langle
G_{0}\right\rangle$=0.8, 1.0, and 1.2 while the variation around $\left\langle
v^{\ast}\right\rangle$ is fixed at ±5% of $\left\langle
v^{\ast}\right\rangle$=1. (The black, red and blue color represents ±5%, ±10%,
and ±20% variation respectively. The dots corresponds to Monte Carlo
simulation results, the line corresponds to the extended model and the dotted
line correspond to the linear fit.)

In Figure 6, the cumulative distributions are actually plotted in form of the
double-logarithm of reciprocal of survival probability, $1-\Phi(\sigma)$,
similar to how it is shown in Figure 2. It must be recalled that this double-
logarithmic representation used in Figure 2 here and numerous other studies
reported by various authors are inspired by the form of Eq. (2) involving a
double-exponential function, which leads to the linear behavior as verified in
Sec. 3. Nevertheless, Eq. (2) itself was derived with the assumption of
constant activation parameters. In contrast, the cumulative distribution,
$\Phi(\sigma)$, given in Eqs. (10) and (12) are observed to be the expectation
of the distribution function, $F(\sigma)$, given in Eq. (2), over the ensemble
of all the possible activation parameters.

Unlike Eq. (2), Eq. (12) gives us no reason to expect a linear behavior with
double-logarithmic scaling and Figs. 6 indeed reflects the non-linearity that
is introduced by allowing the activation parameters to have some variance
around their central values. The scattered plots in Figure 6 illustrate the
results of the Monte-Carlo calculations. We find that by increasing the range
over which the intercept barrier $\left\langle G_{0}\right\rangle$ and
activation volume $v^{\ast}$ are allowed to vary, the extent of curvature of
the plot can be increased. Moreover, cumulative distributions are separately
computed using Eq. (12) numerically, and are also shown with the red solid
curves. It can be seen that the Monte-Carlo results are in excellent agreement
with the analytically derived expression of Eq. (12).

###  IV.4 Implications of parametric variance

The nonlinear trend witnessed in Figure 6 is not limited to our Monte-Carlo
calculations but is often observed in experimental studies as well. In the
original work by Schuh and Lund, from which the statistical model originated,
the authors acknowledged that the tails of the cumulative distribution did not
fit well with the analytical expression, and the linear fit was obtained for
the region excluding the tails (Figs. 5 and 6 in Ref. [21]). Akin to the
results shown in Figure 5, curvatures in the experimentally obtained double-
logarithmic plots have been observed not only in their subsequent study (Fig.
5 in Ref. [22]) but several other experimental and simulation studies have
shown that the tails of the measured distribution never fit to the analytical
solution (_e.g._ , Fig. 1 in Ref. [33]). Comparing Figs. 2 and 6 reveals that
the non-linearity is visible in the present MD simulation results as well.

In all the reported studies, the above-discussed issue is typically bypassed
by selecting a “convenient” part of the whole distribution for fitting. To
see, whether this approach can provide an acceptable estimate of the
activation parameters, we select apparently linear regions in the plots of
Figure. 6 and attempt linear fits (green dashed lines). The resulting
activation parameters computed on the basis of Eq. (3) are listed in Table II,
which reveals an astonishing trend. We find that the activation volumes
$v^{*}_{\text{lin}}$ measured in this manner are much smaller the
corresponding central values, $\left\langle v^{\ast}\right\rangle$, as used in
the Monte-Carlo calculations. Furthermore, the measured values drop even far
below the actual range of $v^{\ast}$. For instance, consider the case of
$\left\langle v^{\ast}\right\rangle$=0.8 and a variation of ±20% around this
mean value, thus allowing $v^{\ast}$ to vary from 0.64 to 0.96. For ±5%
variation in $G_{0}$ around $\left\langle G_{0}\right\rangle$ = 1, the
activation volume obtained using the linear fit is merely 0.2587, which is
only a small fraction of the actual central value, $\left\langle
v^{\ast}\right\rangle$, but also much smaller than the lower bound of 0.85485.
Hence, the naive approach based on Eq. (3) is observed to cause severe
underestimation of the activation volume even for small variances of the
activation parameters. A similar trend is observed in Table III for the
measured values,$G^{0}_{\text{lin}}$ , which are obtained on the basis of Eq.
(3). We find that, similar to the case of actual vs. measured activation
volumes, $G^{0}_{\text{lin}}$ can turn out to be significantly smaller than
$\left\langle G_{0}\right\rangle$ used in the Monte Carlo calculations.

Table 2: $\left\langle G_{0}\right\rangle$ = 1 with $\Delta G_{0}$ as 5% of $\left\langle G_{0}\right\rangle$ and variations of 0%, $\pm$5%, $\pm$10%, and $\pm$20% around $\left\langle v^{\ast}\right\rangle$ $\left\langle v^{\ast}\right\rangle$ | $v^{*}_{\text{lin}}$ |
---|---|---
| 0% | 5% | 10% | 20% |
0.8 | 0.85485 | 0.708 | 0.585 | 0.2587 |
1 | 1.111 | 0.8462 | 0.5965 | 0.3112 |
1.2 | 1.3035 | 1.2165 | 0.6003 | 0.3793 |
Table 3: $\left\langle v^{\ast}\right\rangle$ = 1 with $\Delta v^{\ast}$ as 5% of $\left\langle v^{\ast}\right\rangle$ and variations of 0%, $\pm$5%, $\pm$10%, and $\pm$20% around $\left\langle G_{0}\right\rangle$ $\left\langle G_{0}\right\rangle$ | $G^{0}_{\text{lin}}$ |
---|---|---
| 0% | 5% | 10% | 20% |
0.8 | 0.8828 | 0.8544 | 0.6032 | 0.3676 |
1 | 1.111 | 0.86 | 0.5285 | 0.3806 |
1.2 | 1.0236 | 0.8106 | 0.576 | 0.351 |

The preceding findings suggest that non-trivial mechanisms are not necessarily
necessary to explain the ultra-small activation volumes reported in the
experimental and simulation studies. The statistical method employed in such
studies is mathematically valid only under the assumption of fixed activation
parameters with no variance. As discussed in Sec. IV(B), the activation
parameters are likely to exhibit some variation during repeated measurements
of the critical loads due to a variety of reasons. It is evident from Figure 6
and Tables II and III that the variation of only a few percent is enough to
suppress the measured $v^{\ast}$ far below the actual range of its variation.
It suggests that the estimated activation volumes of the order of atomic
volume are merely the outcome of the premise of the statistical model that is
inapt against the backdrop of physical reality. The actual activation volume
can be much larger owing to the variance in the activation parameters.

##  V Conclusions

In a nutshell, the present work closely examines the strategy of estimating
the activation volume of critical events like dislocation nucleation using
statistical events by analyzing the cumulative distribution of critical loads.
Through the series of measurements, we find that the method remains valid only
in the case of rigidly fixed activation parameters, while the variance in the
measured critical loads is solely attributed to thermal activation. We argue
that this constraint is physically unreasonable, and the activation parameters
are expected to vary over some range. We extend and generalize the statistical
model to allow for such statistical parametric variation and find that even a
little variance in the parameters can significantly alter the distribution of
critical loads in comparison to the conventional constrained model. In
particular, using the method of constrained model on the results obtained with
parametric fluctuations produces the misleading result of extremely small
activation volumes, which explains why the experiments and simulations have
ubiquitously reported activation volumes comparable to atomic volumes. While
the study underlines the effect of not incorporating the parametric
fluctuations within the purview of the constrained statistical model, a major
issue remains unresolved. A generalization of the constrained model to its
unconstrained version shows that the cumulative distribution of critical loads
depends on the probability density of the activation parameters. In the
present study, this is verified by considering a toy distribution, where the
parameters are uniformly distributed over a narrow range. This assumption,
though simple, captures two essential components of a distribution, namely, a
central value and a measure of dispersion. For a realistic physical scenario,
either in experiment or simulation, we do not have any information regarding
the exact form of the parametric distribution. A possible workaround is to
assume a reasonable distribution like a Gaussian or log-normal distribution
and treat its mean and standard deviation as free parameters. Subsequently,
the fitting can be performed using the data from multiple sets of measurements
at various temperatures and loading rates by tuning the distribution
parameters. We hope the outcome of this work can pave the way for further
investigations in this direction.

## References

  * Pöhl [2019] F. Pöhl, Pop-in behavior and elastic-to-plastic transition of polycrystalline pure iron during sharp nanoindentation, Scientific reports 9, 1 (2019).
  * Jiapeng  _et al._ [2017] S. Jiapeng, L. Cheng, J. Han, A. Ma, and L. Fang, Nanoindentation induced deformation and pop-in events in a silicon crystal: molecular dynamics simulation and experiment, Scientific reports 7, 10282 (2017).
  * Ohmura and Wakeda [2021] T. Ohmura and M. Wakeda, Pop-in phenomenon as a fundamental plasticity probed by nanoindentation technique, Materials 14, 1879 (2021).
  * Huang  _et al._ [2009] Y. Huang, Y. L. Chiu, J. Shen, J. J. Chen, and J. Sun, Nanoindentation study of ti-based metallic glasses, Journal of alloys and compounds 479, 121 (2009).
  * Mojumder  _et al._ [2020] S. Mojumder, M. Mahboob, and M. Motalab, Atomistic and finite element study of nanoindentation in pure aluminum, Materials Today Communications 23, 100798 (2020).
  * Zhao  _et al._ [2018] K. Zhao, A. Mayer, J. He, and Z. Zhang, Dislocation based plasticity in the case of nanoindentation, International Journal of Mechanical Sciences 148, 158 (2018).
  * Lee  _et al._ [2020] S. Lee, A. Vaid, J. Im, B. Kim, A. Prakash, J. Guénolé, D. Kiener, E. Bitzek, and S. H. Oh, In-situ observation of the initiation of plasticity by nucleation of prismatic dislocation loops, Nature communications 11, 2367 (2020).
  * Dahmen  _et al._ [2011] K. A. Dahmen, Y. Ben-Zion, and J. T. Uhl, A simple analytic theory for the statistics of avalanches in sheared granular materials, Nature Physics 7, 554 (2011).
  * Friedman  _et al._ [2012] N. Friedman, A. T. Jennings, G. Tsekenis, J.-Y. Kim, M. Tao, J. T. Uhl, J. R. Greer, and K. A. Dahmen, Statistics of dislocation slip avalanches in nanosized single crystals show tuned critical behavior predicted by a simple mean field model, Physical review letters 109, 095507 (2012).
  * Avila  _et al._ [2019] K. E. Avila, S. Küchemann, I. Alabd Alhafez, and H. M. Urbassek, Shear-transformation zone activation during loading and unloading in nanoindentation of metallic glasses, Materials 12, 1477 (2019).
  * Yu  _et al._ [2021] J. Yu, A. Datye, Z. Chen, C. Zhou, O. E. Dagdeviren, J. Schroers, and U. D. Schwarz, Atomic-scale homogeneous plastic flow beyond near-theoretical yield stress in a metallic glass, Communications Materials 2, 22 (2021).
  * Rao  _et al._ [2011] S. Rao, D. Dimiduk, T. Parthasarathy, J. El-Awady, C. Woodward, and M. Uchic, Calculations of intersection cross-slip activation energies in fcc metals using nudged elastic band method, Acta materialia 59, 7135 (2011).
  * Geslin  _et al._ [2017] P.-A. Geslin, R. Gatti, B. Devincre, and D. Rodney, Implementation of the nudged elastic band method in a dislocation dynamics formalism: Application to dislocation nucleation, Journal of the Mechanics and Physics of Solids 108, 49 (2017).
  * McPhie  _et al._ [2012] M. McPhie, S. Berbenni, and M. Cherkaoui, Activation energy for nucleation of partial dislocation from grain boundaries, Computational materials science 62, 169 (2012).
  * Ryu  _et al._ [2011a] S. Ryu, K. Kang, and W. Cai, Predicting the dislocation nucleation rate as a function of temperature and stress, Journal of Materials Research 26, 2335 (2011a).
  * Ryu  _et al._ [2011b] S. Ryu, K. Kang, and W. Cai, Entropic effect on the rate of dislocation nucleation, Proceedings of the National Academy of Sciences 108, 5174 (2011b).
  * Weinan  _et al._ [2005] E. Weinan, W. Ren, E. Vanden-Eijnden, _et al._ , Finite temperature string method for the study of rare events, J. Phys. Chem. B 109, 6688 (2005).
  * Aziz  _et al._ [2006] M. J. Aziz, Y. Zhao, H.-J. Gossmann, S. Mitha, S. P. Smith, and D. Schiferl, Pressure and stress effects on the diffusion of b and sb in si and si-ge alloys, Physical Review B 73, 054101 (2006).
  * Sahu  _et al._ [2019] B. P. Sahu, A. Dutta, and R. Mitra, Mechanism of negative strain rate sensitivity in metallic glass film, Journal of Alloys and Compounds 784, 488 (2019).
  * Kato [2009] M. Kato, Thermally activated dislocation depinning at a grain boundary in nanocrystalline and ultrafine-grained materials, Materials Science and Engineering: A 516, 276 (2009).
  * Schuh and Lund [2004] C. A. Schuh and A. C. Lund, Application of nucleation theory to the rate dependence of incipient plasticity during nanoindentation, Journal of Materials research 19, 2152 (2004).
  * Mason  _et al._ [2006] J. Mason, A. Lund, and C. Schuh, Determining the activation energy and volume for the onset of plasticity during nanoindentation, Physical review B 73, 054102 (2006).
  * Mendelev  _et al._ [2003] M. Mendelev, S. Han, D. Srolovitz, G. Ackland, D. Sun, and M. Asta, Development of new interatomic potentials appropriate for crystalline and liquid iron, Philosophical magazine 83, 3977 (2003).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, _et al._ , Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271, 108171 (2022).
  * Stukowski [2009] A. Stukowski, Visualization and analysis of atomistic simulation data with ovito–the open visualization tool, Modelling and simulation in materials science and engineering 18, 015012 (2009).
  * Dutta [2017] A. Dutta, Compressive deformation of fe nanopillar at high strain rate: modalities of dislocation dynamics, Acta Materialia 125, 219 (2017).
  * Zhu  _et al._ [2013] C. Zhu, Z. Lu, and T. Nieh, Incipient plasticity and dislocation nucleation of fecocrnimn high-entropy alloy, Acta Materialia 61, 2993 (2013).
  * Ye  _et al._ [2017] Y. Ye, Z. Lu, and T. Nieh, Dislocation nucleation during nanoindentation in a body-centered cubic tizrhfnb high-entropy alloy, Scripta Materialia 130, 64 (2017).
  * Chachamovitz and Mordehai [2018] D. Chachamovitz and D. Mordehai, The stress-dependent activation parameters for dislocation nucleation in molybdenum nanoparticles, Scientific reports 8, 3915 (2018).
  * Brandl  _et al._ [2022] C. Brandl, I.-C. Choi, and R. Schwaiger, Temperature-dependent nanoindentation and activation volume in high-purity body-centered cubic chromium, Materials Science and Engineering: A 852, 143629 (2022).
  * Tan  _et al._ [2023] A. M. Z. Tan, Z. Li, and H. Gao, Statistical modeling of the effect of chemical inhomogeneity on incipient plasticity in complex concentrated alloys, Physical Review Materials 7, 053601 (2023).
  * Nisany and Mordehai [2022] S. Nisany and D. Mordehai, A multiple site type nucleation model and its application to the probabilistic strength of pd nanowires, Metals 12, 280 (2022).
  * Derlet and Maass [2016] P. Derlet and R. Maass, The stress statistics of the first pop-in or discrete plastic event in crystal plasticity, Journal of Applied Physics 120, 225101 (2016).