# A stochastic discrete slip approach to microplasticity: Application to
submicron W pillars

Carlos J. Ruestes  Javier Segurado

###### Abstract

A stochastic discrete slip approach is proposed to model plastic deformation
in submicron domains. The model is applied to the study of submicron pillar
($D~{}\leq~{}1\mu m$) compression experiments on tungsten (W), a prototypical
metal for applications under extreme conditions. Slip events are geometrically
resolved in the specimen and considered as eigenstrain fields producing a
displacement jump across a slip plane. This novel method includes several
aspects of utmost importance to small-scale plasticity, i.e. source truncation
effects, surface nucleation effects, starvation effects, slip localization and
an inherently stochastic response. Implementation on an FFT-spectral solver
results in an efficient computational 3-D framework. Simulations of submicron
W pillars ($D~{}\leq~{}1\mu m$) under compression show that the method is
capable of capturing salient features of sub-micron scale plasticity. These
include the natural competition between pre-existing dislocations and surface
nucleation of new dislocations. Our results predict distinctive flow stress
power-law dependence exponents as well as a size-dependence of the strain-rate
sensitivity exponent. The results are thoroughly compared with experimental
literature.

###### keywords:

plasticity , simulation , slip localization , stochasticity , tungsten

††journal: International Journal of Plasticity

[inst1]organization=IMDEA Materials Institute,addressline=C/Eric Kandel 2,
city=Getafe, postcode=28906, state=Madrid, country=Spain

[inst2]organization=Instituto Interdisciplinario de Ciencias Básicas (ICB),
Universidad Nacional de Cuyo UNCuyo-CONICET, Facultad de Ciencias Exactas y
Naturales,addressline=Padre Contreras 1300, city=Mendoza, postcode=5500,
state=Mendoza, country=Argentina

[inst3]organization=Universidad Politécnica de Madrid, Department of Materials
Science,addressline=E.T.S.I. Caminos, city=Madrid, postcode=28040,
state=Madrid, country=Spain

##  1 Introduction

Progress towards fulfilling the promise of nanoscale solutions to engineering
problems hinges on the development of computational and characterization
techniques for an adequate assessment of the mechanical response of the
resulting materials and parts. For instance, micropillar compression testing
(Uchic et al., 2004; Uchic and Dimiduk, 2005) can provide critical data for
the development of a number of technologies at different engineering scales,
from the micro/nanoscale (e.g micro/nano electromechanical systems -MEMS /
NEMS-, microsensors and actuators and microscale energy harvesting and storage
systems) to the macroscale. The technique is instrumental in establishing
correlations between the microstructure and the mechanical properties of
materials (Greer and De Hosson, 2011; Dehm et al., 2018). The analysis and
interpretation of the experimental data obtained often relies on advanced
computational simulations at different length scales, which aim at determining
the underlying deformation micromechanisms responsible of the observed
behavior.

At the lowest scale, molecular dynamics simulations (MD) were used to study
the interaction of dislocations with pillar surfaces, leading to the proposal
of a surface-induced cross-slip mechanism (Weinberger and Cai, 2008). In
addition, several studies have been published with different focuses, such as
unraveling the tension-compression asymmetry observed experimentally (Healy
and Ackland, 2014) or analyzing the origin of size effects observed
experimentally (Yang et al., 2021; Yaghoobi and Voyiadjis, 2016), to name a
few. Nevertheless, the detailed nature of atomistic simulations, together with
integration timesteps of the order of 1-3 fs, enforces the use of ultra-high
strain rates of $10^{7}/s$ and above and very small pillar sizes
($D~{}<~{}100~{}nm$). The comparison between MD simulations and micropillar
compression experiments is thus limited by the inherent drawbacks associated
with the simulation technique. Still, MD has made fundamental contributions to
the understanding of source nucleation in micropillars and nanowires
(Weinberger and Cai, 2008; Zhu et al., 2008), even allowing to estimate
nucleation stresses on pristine fcc nanowires within 20% of the experimental
values (Jennings et al., 2013).

At the mesoscale, discrete dislocation dynamics simulations (DDD) have been
used to provide insights into the mechanisms responsible for the
experimentally observed staircase stress-strain behavior (Tang et al., 2008)
as well as athermal size-dependent strengthening (Rao et al., 2008). The
technique has also been used to study the effect of specimen size on plastic
flow and work-hardening (El-Awady et al., 2009) and also to study the
characteristics of source truncation controlled flow behavior for submicron
FCC single crystal micropillars (Cui et al., 2014). Moreover, in order to
capture surface-related effects, MD simulation-derived models of surface
nucleation (Ryu et al., 2011; Zhu et al., 2008) and surface-induced cross-slip
(Weinberger and Cai, 2008) have been incorporated into DDD frameworks (Ryu et
al., 2013, 2015; Hu et al., 2019), allowing for the exploration of pillar
sizes in the range of 100 nm to a few microns. The technique has also been
used to develop stochastic models for the onset of plasticity in micro- and
nano-scale structures (Shao et al., 2014). In spite of their advantages, DDD
simulations still suffer from timescale limitations (few ns to few $\mu$s) and
typical strain rates used are in the range of $10^{3}/s-10^{6}/s$, though
strain rates as small as 0.1/s are possible depending on the computational
resources (Fan et al., 2021). Both techniques, MD and DDD, typically rely on
an implementation in the form of massively parallel solvers to obtain results
in a reasonable time, which in turn impose important hardware requirements.

At the continuum scale, several groups had opted for finite element
simulations (FE) with different material model approaches. FE simulations
using simple isotropic plasticity have proven useful in providing
recommendations for experiments design and reliable testing (Zhang et al.,
2006; Kiener et al., 2009). Simulations using crystal plasticity constitutive
equations (CPFE) simulations have been successfully used to study the effects
of initial crystallographic orientation, diameter-to-length ratio and friction
on the response of Cu single crystalline pillars (Raabe et al., 2007). CPFE
has also been used for the study of polycrystalline Ni-based superalloy
micropillars (Cruzado et al., 2015). The use of CPFE models has also extended
to size effects in Al pillar compression, focusing on a continuum description
of starvation (Jérusalem et al., 2012). CPFE relies on homogenization of the
evolution of large dislocation ensembles and on scale-separation, hypotheses
that can be perfectly assumed when studying the deformation of large crystal
grains or specimens of several microns. However, in small-scale testing the
number of dislocations in the specimen can be very small (just a few
dislocations in the full specimen), such that their effect on the deformation
is localized in specific planes and show an stochastic nature. These aspects
cannot be homogenized in a continuum-scale law.

In order to fill gaps between DDD and CPFE models several attempts have been
made to extend standard CPFE models to account for some characteristic
features of small-size tests (i.e. the stochastic nature of the response and
slip-band localization). Regarding the stochastic response, Ng and Ngan (2008)
proposed a Monte Carlo model based on the survival probability and the burst
size vs. stress distributions but with no further connection with the material
structure. Later on, Konstantinidis et al. (2014) proposed a cellular
automaton, based on the gradient plasticity framework of Zhang and Aifantis
(2011), to model stochastic effects in pillar compression. In addition, Lin et
al. (2015) presented a crystal plasticity model in which stochasticity was
introduced by defining the plastic strain as composed of a series of strain
bursts whose sizes follow a power-law distribution function and whose rates
are determined by a constitutive equation. More recently, a general stochastic
CP formulation based on kinetic Monte Carlo, was proposed in which the active
slip systems were selected based on the strain rates dictated by a mobility
law depending on resolved stress and temperature (Yu et al., 2021). However,
the homogenized nature of the formulation prevents its direct use in
applications where slip localization is expected.

The slip-band localization observed in experiments, resulting from the
activation of discrete slip events in a particular plane, cannot be naturally
captured by CPFE models. Some attempts have been made to emulate this
localization in CP frameworks. Lin et al. (2016) developed a dislocation-based
crystal plasticity model in which, in order to trigger localization in a
single slip band, strain-softening was added to the slip resistance evolution
law and a void was introduced in the center of the pillar to trigger
localization. More recently, Wijnen et al. (2021) presented a CPFE model
adapted for simulating heterogeneous plastic deformation in single crystals.
In this model, only some regions of the domain under study are modeled using
CP keeping the others as elastic material. These slip areas are sampled from a
probability distribution based on characteristics of the dislocation sources
and plastic localization is accomplished by considering different slip bands
and the weakest link principle (Norfleet et al., 2008). These extended CPFE
models allow to describe the strain localization in pillars but maintain a
deterministic nature and do not explicitly consider plastic events by
individual dislocations.

Summarizing, the mechanical response of metals and alloys under micro
compression testing is typically characterized by certain features: i) samples
usually have a low dislocation density and exhibit source truncation
hardening; ii) Plastic deformation is usually localized and in the form of
slip bands; iii) For diameters below 150 nm, or starved samples such as those
produced by mechanical annealing, surface nucleation of dislocations is likely
to occur. As a consequence, plasticity takes place by events of discrete
nature and the mechanical response of such structures is typically stochastic.
Despite all the aforementioned computational advances, capturing all these
features into a single and efficient computationally-cheap 3-D framework is
still a challenge.

The purpose of this work is two fold. Firstly, we present a multiscale
computational framework capable of capturing the fundamental aspects of the
mechanical response of submicron metallic parts. In this model, the
displacement jumps along a cross sectional plane of the sample produced by the
slip of a dislocation are introduced as eigenstrains around that plane using
the Eshelby formalism. The framework also accounts for surface nucleation by
implementing a model derived from MD simulations. The stochastic nature of the
mechanical response is captured by means of a kinetic Monte Carlo selection
process. In this kMC approach, the events correspond to the single slip events
produced by dislocations, whose probabilities are dictated by physically-based
laws. This approach is implemented in an FFT-based solver resulting in a very
efficient and computationally-cheap framework that allows to model complex 3D
geometries. Secondly, we revisit micropillar compression experiments of single
crystalline bcc submicron pillars ($D~{}\leq~{}1\mu m$) with focus on tungsten
(W). We show that not only our method can successfully capture fundamental
aspects of the mechanical response at these scales, but also allows for
detailed investigation of size-effects, strain-rate effects and statistics of
strain-bursts.

##  2 Methodology

The framework developed assumes dislocation-mediated plasticity as the
underlying inelastic deformation micromechanism. We consider pre-existing
dislocations as well as newly developed dislocations due to surface
nucleation. The amount of pre-existing dislocations $N_{d}$ on a pillar of
diameter $D=2R$, height $H$, and initial dislocation density $\rho$, is
determined as

$N_{d}=int(\rho~{}\frac{\pi~{}R^{2}~{}H}{L_{d}})$ (1)

where $L_{d}$ is the average length of a dislocation segment, here taken as R
for simplicity. Then, $N_{d}$ slip planes are randomly selected along the
pillar axis with normal directions consistent with the available slip systems.

Surface nucleation can potentially happen at any position on the pillar
surface. From a practical point of view, in order to consider a dislocation
nucleation and its subsequent slip in our framework, a large but finite number
of nucleation sites $N_{SN}$ is considered. In practice, $100<N_{SN}<200$. The
positions of the resulting slip planes are randomly selected along the pillar
axis, and their orientations are also based on the available slip systems.

In our framework, dislocations are not explicitly considered, but the effect
of their extension when reaching the surface of the pillar. This effect
consists in a relative displacement between the upper and lower parts of their
plane and experimentally would define a slip trace. This description is
equivalent to the definition of a closed dislocation loop in phase-field
dislocation dynamics as an inclusion with an eigenstrain (an Eshelby
inclusion) (Wang et al., 2001; Rodney et al., 2003). In our case, only the
final stage of the loop is considered which consists then in a loop with a
shape defined by the intersection of the pillar with the slip plane; see Sec.
2.1 for details. The different parts of the model are introduced in the next
sections.

###  2.1 Kinematics of a slip event

First we derive the elastic displacement introduced by a dislocation loop
considering it as an equivalent Eshelbian inclusion. Under this
representation, the loop is idealized by moving the upper side of surface $S$,
denoted by $S^{+}$, by $b$ with respect to the lower side $S^{-}$, as shown in
Fig. 1.a. Let $\mathbf{n}$ be the normal to the loop, then the relative
displacement between the upper and lower parts of the loop in direction
$\mathbf{n}$ is given by

$\llbracket\mathbf{u}\rrbracket_{n}(\mathbf{x})=\mathbf{b}\
\delta_{S}(\mathbf{x}),$ (2)

where $\mathbf{b}$ is the burgers vector of the loop and
$\delta_{S}(\mathbf{x})$ denotes the surface Dirac delta function, that is
infinite if $\mathbf{x}\in S$ and zero elsewhere. If the loop is fully
embedded in a material, the resulting elastic fields of this incompatible
deformation can be obtained using Green’s function (Mura, 2013). In our model,
the only stage considered is when the loop has filled all the pillar cross
section.

The displacement jump due to a full plane slip has been considered in the past
using a strong discontinuity approach in the context of finite elements and
discrete dislocation dynamics (Romero et al., 2008) . However, introducing
this discontinuity in other numerical frameworks as FFT imply smoothing out
the jump to make the model numerically tractable. This can be done by
distributing the jump along a thin band of size $h$, such that it renders the
same total displacement as in the initial ideal case. In a simplified way, the
line integral along a path parallel to $\mathbf{n}$ corresponds to

$\int_{z=-h/2}^{z=h/2}\boldsymbol{\varepsilon}(\mathbf{x})\mathrm{d}z=\llbracket\mathbf{u}\rrbracket_{n}$
(3)

Figure 1: Dislocation eigenstrain formalism and the application to a single
crystalline pillar. (a) A dislocation loop $L$ acting as a boundary of a
dislocation surface $S$. By sweeping the upper part $S^{+}$ with respect to
$S^{-}$ in an amount $b$, a dislocation is effectively introduced. (b) For a
dislocation sweeping a pillar, the dislocation can be replaced by a plate-like
Eshelbian inclusion of thickness h.

Extending the use of the equivalent Eshelbian inclusion, Figure 1.b depicts a
scenario in which the cross-section of cylindrical pillar has been fully swept
in an amount of $b$ by a single dislocation loop. The associated displacement
is linearly distributed over a thickness $h$ in the direction perpendicular to
the loop. As explained later on (Sec. 2.4), the simulation domain is
discretized in a box of [Nx, Ny, Nz] voxels, so in practice, $h$ is taken as
twice the distance between voxels (Capolungo and Taupin, 2019). As presented
in the Supplementary Material file, the method is relatively insensitive to
the election of the voxel size. In consequence, the strain
$\boldsymbol{\varepsilon}^{EIG~{}(i)}(\mathbf{x})$ associated with a slip
event $i$, is taken as:

$\boldsymbol{\varepsilon}^{EIG~{}(i)}(\mathbf{x})=\begin{cases}\frac{b}{h}~{}(\boldsymbol{s}^{i}\otimes\boldsymbol{n}^{i}),~{}if~{}\mathbf{x}~{}\epsilon~{}plastic~{}region\\\
0,~{}elsewhere\end{cases}$ (4)

where $b$ is the burgers vector modulus, $h$ is the thickness of the slip band
and $\mathbf{s}^{i}$ and $\mathbf{n}^{i}$ stand for the burgers vector
direction (slip) and normal direction to the slip plane,
$\mathbf{s}^{i}\otimes\mathbf{n}^{i}$ being the Schmid tensor. Consequently,
for a pillar that underwent $M$ slip events, the total eigenstrain field
$\boldsymbol{\varepsilon}^{EIG}(\mathbf{x})$ is the superposition of all the
eigenstrains applied,

$\boldsymbol{\varepsilon}^{EIG}(\mathbf{x})=\sum_{i=1}^{M}\boldsymbol{\varepsilon}^{EIG~{}(i)}(\mathbf{x})$
(5)

###  2.2 Constitutive laws for the slip events

The laws linking the stress at a given plane with its slip rate will be
presented in next sections.

####  2.2.1 Pre-existing dislocations

The model assumes that the pillar contains a certain number of dislocations
which will drive the slip events. Two types of pre-existing dislocations are
considered, these are single-arm sources and screw dislocation segments. The
former aid in the introduction of source truncation effects, while the latter
allow for source exhaustion effects.

##### Single-Arm Sources

A dislocation with a pinning point in the interior of the crystal and a free
end on the surface acts as a source of dislocations in small pillars and is
commonly termed as Single-Arm Source (SAS) (Parthasarathy et al., 2007). In
the present model, these sources are randomly distributed in the pillar
volume, assigning to each source a particular slip plane. The critical
resolved stress for activating a SAS is determined by its arm length
$\lambda$,

$\tau^{SAS}_{CRSS}=\frac{\alpha_{SAS}~{}\mu~{}b}{\lambda}+\tau_{0}$ (6)

where $\alpha_{SAS}$ is a source-strength coefficient that considers the
nature of the source, its length and material, here taken as 0.6 after (Rao et
al., 2007). $\mu$ is the shear modulus, $b$ the Burgers vector and $\tau_{0}$
the friction stress. In their presentation of the SAS model, Parthasarathy et
al. (2007) include a Taylor-type term in eq. 6 to account for dislocation
forest interactions. This term is not considered here since, for D $\leq$ 1
$\mu$m specimens and typical dislocation densities of the order of
$10^{12}m^{-2}-10^{13}m^{-2}$, the number of dislocations present in the
pillar becomes too small to homogenize their effect ($N_{d}~{}\leq~{}10)$.
Under such conditions, weak interactions are expected. Thus, we essentially
neglect hardening. For pillar diameters above 1 micron, hardening should be
incorporated, as explained in Sec. 4.

In order to generate a random arrangement of SAS, for each source a random
pining point is generated in the pillar interior. Then, a slip plane among the
characteristic active planes of the metal lattice is randomly assigned. The
source is defined as a straight line of length $\lambda$ connecting the
pinning point with the nearest point in the elliptical boundary resulting from
the intersection of the plane with the pillar surface. For more details on
this model, the reader is referred to (Parthasarathy et al., 2007).

This model has been widely used on experimental studies (Abad et al., 2016),
DDD simulations (Hu et al., 2019), Crystal Plasticity simulations (Gu et al.,
2021), as well as on recent discrete slip plane models (Wijnen et al., 2021).
In the present work, it is assumed that single-arm sources remain immutable
after activation. Changes in SAS could be potentially included in the model,
for example the incorporation of SAS destruction as proposed in the work by
Cui et al. (2015) on coated submicron pillars.

##### Screw dislocation segments

In addition to SASs, we also consider the possibility of pre-existing pure-
screw segments. In essence, such dislocations are not pinned and can sweep
across the slip plane, exiting the cylinder as they reach its surface. Under
these operating conditions, screw segments allow to capture the so-called
”dislocation-starvation” behavior (Greer and Nix, 2006). In this case, the
resolved stress necessary to move the dislocation is only the lattice
friction,

$\tau^{screw}_{CRSS}=\tau_{0}$ (7)

The position of the screw dislocation within the pillar cross-section is
randomly chosen. As shown in section 3.3, this allows capturing different
cases of mechanical annealing observed in the literature, with variable stress
drop due to the activation of these sources. It is worth noting that in our
implementation, dislocation cross-slip due to image forces in the vicinity of
the free surfaces is not considered. This is further discussed in Sec. 4.

####  2.2.2 Surface nucleation model

For pillars and nanowires in the range of a few tenth of nanometers to a few
hundred of nanometers, virtually defect-free microstructures are not unlikely.
Under these conditions, mechanical testing of nanopillars and nanowires reveal
extremely high tensile stress and atomistic simulations point to surface
nucleation of dislocations as the source of plasticity (Zhu et al., 2008). In
consequence, the framework developed includes surface nucleation effects.

The surface nucleation model adopted here was initially proposed after
atomistic simulations (Zhu et al., 2008; Ryu et al., 2011) and later adopted
on discrete dislocation dynamics (DDD) simulations (Ryu et al., 2015; Hu et
al., 2017, 2019). Briefly, within a time span $\Delta t$, the probability of a
surface nucleation event in a pillar of surface $S$ is given by

$P=\nu_{0}~{}exp[-\frac{F(\sigma,T)}{k_{B}~{}T}]\cdot\frac{S}{b^{2}}\cdot\Delta
t.$ (8)

Here, $\nu_{0}$ is an attempt frequency and $k_{B}$ the Boltzmann
constant.Thus, the ratio $\frac{S}{b^{2}}$ stands for the number of potential
nucleation sites. Finally, $F(\sigma,T)$ is an activation free energy that
depends on stress and temperature through

$F(\sigma,T)=(1-\frac{T}{T_{d}})\cdot F_{0}(\sigma)$ (9)

with

$F_{0}(\sigma)=A\cdot(1-\frac{\sigma}{\sigma_{athm}})^{\alpha_{SN}}.$ (10)

Here, $A$ is the zero-stress, zero-temperature activation energy, while
$\alpha_{SN}$ and $T_{d}$ are constants corresponding to the stress and
temperature dependencies, respectively. $T_{d}$ is usually taken as half the
bulk melting temperature (Zhu et al., 2008). Finally, $\sigma_{athm}$ is an
athermal nucleation stress, corresponding to the tensile stress required for
surface nucleation at zero-temperature. The possibility of surface nucleation
is considered once $P=1$.

The parameters $A$, $T_{d}$, $\alpha_{SN}$ and $\sigma_{athm}$ are material-
dependent and their calibration require dedicated studies, either based on
atomistic simulations (Zhu et al., 2008) or on experiments (Chen et al.,
2015). To the best of the authors’ knowledge, there are no such studies on W.
Therefore, we opt to fit these parameters to the available experimental data,
following suggestions by Chen et al. (2015). Nanopillar compression and
nanowire tensile test for [100]-oriented W are extremely scarce in the
literature. Srivastava et al. (2021) report yield stresses of the order of 3.5
GPa in their 100 nm diameter pillar compression tests, whereas Córdoba et al.
(2017) report yield stresses of 6.6 $\pm$ 1.2 GPa for their 93 nm suspended
wires under bending. Therefore, we opted to fit the parameters of eqn. 8 to a
100 nm target yield strength of 4.5 GPa, in between the mean value of
Srivastava et al. (2021) and the lower bound of the range reported by Córdoba
et al. (2017). The resulting parameters are included in Table 1 after the
fitting procedure presented in the Supplementary Material file. The surface
nucleated dislocations are considered as pure screw segments, as described in
the previous subsection.

Symbol | Property | Value | Unit
---|---|---|---
$T_{d}$ | Disorder temperature | 1800 | K
A | zero $\sigma$, zero $T$ activation energy | 6 | eV
$\alpha_{SN}$ | Fitting exponent stress dependence | 4 | -
$\sigma_{athm}$ | Athermal nucleation stress | 16 | GPa
Table 1: Surface nucleation model parameters.

####  2.2.3 Dislocation mobility and displacement rates

The displacement rate of the pillar in the slip direction $\mathbf{s}^{i}$ due
to the slip produced in a plane $i$ by the movement of a dislocation can be
obtained using a reasoning similar to Orowan’s equation,

$\dot{\mathbf{u}}^{i}_{p}=b^{i}\frac{v^{i}}{D^{i}}\mathbf{s}^{i}$ (11)

where $b^{i}$ and $v^{i}$ are the Burgers vector modulus, and dislocation
velocity corresponding to source $i$, respectively. $D^{i}$ corresponds to the
length in direction $\mathbf{s}^{i}$ swept by a single dislocation when fully
crossing the slip plane $i$. If the unit vector $\mathbf{e}$ stands for the
pillar axis and $D$ for its diameter, then
$D^{i}=D/(\mathbf{n}^{i}\cdot\mathbf{e}$), and the vertical displacement rate
corresponds to

$\dot{u}^{i}_{p}=\dot{\mathbf{u}}^{i}_{p}\cdot\mathbf{e}=b^{i}\frac{v^{i}}{D/\mathbf{n}^{i}\cdot\mathbf{e}}\mathbf{s}^{i}\cdot\mathbf{e}=\frac{b^{i}\
v^{i}M^{i}}{D}$ (12)

where $M^{i}$ is the Schmid factor of slip system associated to source $i$
projected in direction $\mathbf{e}\otimes\mathbf{e}$.

In order to compute the displacement rate for each plane $i$, the associated
dislocation velocity $v^{i}$ must be determined. Dislocation-mediated
plasticity in bcc metals can take place in 48 bcc slip systems including the
$\\{110\\}$,$\\{112\\}$ and $\\{123\\}$ families of slip planes.
Interestingly, some works propose the decomposition of $\\{112\\}$ and
$\\{123\\}$ slip planes on alternating $\\{110\\}$ slip planes (Christian,
1983; Marichal et al., 2013; Ruestes et al., 2014). Based on this, and for the
sake of simplicity, we consider bcc metals with only $\\{110\\}$ slip planes
($\bar{b}=1/2~{}a_{0}~{}<111>$ dislocations). That is, $|\bar{b}|$ is unique
and equal to $a_{0}~{}\sqrt{3}/2$, with $a_{0}$ the bcc lattice parameter.

The mobility laws for bcc tungsten dislocations have been deeply studied using
atomistic simulations (Cereceda et al., 2013; Stukowski et al., 2015; Cereceda
et al., 2016; Po et al., 2016). For tungsten, and considering relatively low
temperatures, dislocation velocities can be estimated by:

$v=\begin{cases}\frac{\mu_{0}~{}h^{\prime}~{}(\zeta-w)}{b}exp(-\frac{\Delta
H(\tau_{eff})}{k_{B}~{}T}),&\Delta H>0\\\ \frac{\tau_{eff}~{}b}{B},&\Delta
H\leq 0\end{cases}$ (13)

$\Delta H(\tau_{eff})=\Delta H_{0}\lceil
1-(\frac{\tau_{eff}}{\tau_{*}})^{p}\rceil^{q}$ (14)

Symbol | Property | Value | Unit
---|---|---|---
$\mu_{0}$ | Attempt frequency | 9.1 $10^{11}$ | $s^{-1}$
$h^{\prime}$ | Distance between Peierls valleys | $\sqrt{6}/3$ | $a_{0}$
$w$ | Kink pair width | 11 | b
$\zeta$ | Mean dislocation segment width | 25 | b
$\tau_{P}$ | Peierls stress | 2.03 | GPa
$\Delta H_{0}$ | Kink-pair energy at 0 K | 1.63 | eV
$p$ | Energy profile parameter | 0.86 | -
$q$ | Energy profile parameter | 1.69 | -
$B$ | Phonon drag coefficient | 9.8 $10^{-4}$ | Pa*s
Table 2: Parameters for the dislocation mobility law. Obtained from molecular
dynamics simulations (Cereceda et al., 2013; Stukowski et al., 2015; Cereceda
et al., 2016; Po et al., 2016)

In the expression, $\mu_{0}$ is an attempt frequency, $h^{\prime}$ is the
distance between Peierls valleys, $\zeta$ is the mean dislocation segment
length for the slip system, $w$ is the kink-pair width and $\Delta{H_{0}}$ is
the kink-pair formation energy. Parameters $p$ and $q$ allow for the
adjustment of the ”tail” ($\tau_{eff}\simeq 0$) and ”top”
($\tau_{eff}\simeq\tau_{*}$) of the energy profile $\Delta H$, respectively
(Kocks et al., 1975). The parameters for W are presented in Table 2.

The upper part of the definition ($\Delta H>0$) corresponds to the kink-pair
thermally-activated regime, typical of bcc metals (Kocks et al., 1975). The
lower part, in turn, corresponds to the phonon-drag-dominated regime, typical
of high stress / high strain-rate conditions (Po et al., 2016). In Eq. (14)
$\tau_{eff}$ represents the effective resolved stress on the plane $i$,

$\tau_{eff}=|\overline{\tau}^{i}|-\tau_{a}^{i}$ (15)

where $\overline{\tau}^{i}$ is the resolved shear stress, averaged over the
plane, for source $i$ resulting from the Cauchy stress projected in the
particular system to which that source belongs to,

$\tau^{i}=\boldsymbol{\sigma}:(\mathbf{s}^{i}\otimes\mathbf{n}^{i}).$ (16)

In Eq. (15) $\tau_{a}$ is the resistance to slip due to athermal barriers, as
dislocation groups, precipitates, which are not considered in the present
study ($\tau_{a}^{i}=0$). Finally, $\tau_{*}$ represents the part of the
resistance to slip due to thermally activatable obstacles. In the case of pure
W the only thermal barrier is the Peierls resistance, hence
$\tau_{*}=\tau_{P}$.

In essence, the use of eq. 13 for the calculation of the pillar displacement
rates (eq. 12) allows to incorporate the thermally-activated motion character
of screw dislocations at low stresses and temperatures, as well as the phonon-
drag-dominated regime for newly surface-nucleated dislocations under high
stresses. Equations 12-13 allow for the effect of each plane $i$ in the
displacement rate of the pillar.

###  2.3 Kinetic Monte Carlo procedure for stochastic slip events

In our model, the activation of the previously described microscopic slip
mechanisms is controlled by a kinetic Monte Carlo selection process.
Considering a displacement-controlled micropillar compression experiment with
an applied strain rate of $\dot{\varepsilon}_{0}$, the velocity of the upper
part of the pillar corresponds to $\dot{u}_{0}=H\dot{\varepsilon}_{0}$. This
displacement is accommodated by either elastic or plastic deformation. The
plastic contribution to the displacement of the pillar corresponds to the
superposition of the displacement produced by slip in the active planes

$\dot{u}_{p}=\sum_{i}\dot{u}^{i}_{p}$

with $\dot{u}^{i}_{p}$ given by eq. 12. The elastic displacement rate then
corresponds to

$\dot{u}_{E}=\dot{u}_{0}-\dot{u}_{p}.$ (17)

Yu et al. (2021) proposed a rejection-free kinetic Monte Carlo algorithm to
obtain the stress-strain evolution in a deforming single crystal in which the
strain rates of the different slip systems provide the set of event rates. In
this work, this idea is generalized to a full deforming specimen solving the
macroscopic stress-strain response as well as the spatial distribution of the
fields involved. At time $t^{n}$, each of the $i$ escape pathways has a rate
constant $\dot{u}_{p}^{i}$, that characterizes the probability per unit time
that the system escapes to that state $i$. These rates can be used to assemble
an array of partial sums representing the accumulated rate of all the objects
up to and including object j,

$r_{j}=\sum_{i=1}^{j}\dot{u}_{p}^{i}+r_{0}$ (18)

with $r_{0}=\dot{\varepsilon_{0}}$ Then, the total escape rate is the sum of
all the rates

$r_{tot}=\sum_{i=1}^{N}\dot{u}_{p}^{i}+r_{0}$ (19)

$N$ being the total number of sliding planes considered. In agreement with
rejection-free kMC theory, the next event to be executed is the $k^{th}$
process satisfying

$r_{k-1}<\xi_{1}r_{tot}\leq r_{k}$ (20)

If the event selected corresponds to the rate $r_{0}$, the next event will
correspond to an elastic event, i.e. no inelastic deformation will be
introduced in the plastic layers. Otherwise, a plastic event is chosen, for
which the corresponding eigenstrain will be applied. The pillar geometry,
array of displacement rates, and the resulting pillar shape after the whole
process are exemplified in Figure 2.

Figure 2: Schematic of stochastic selection process. a) a number of plastic
regions are randomly selected both for pre-existing dislocations (blue) and
possible surface nucleation sites (red). b) Upon loading and based on eqns. 12
through 19, a sampling array of events in built. The total rate $r_{tot}$ is
composed of the prescribed applied strain rate $r_{0}$ and the vertical
displacement rates of each plastic region. In general, the value of each
vertical displacement rate $\dot{u}^{i}_{p}$ is different from the rest of the
displacement rates during each time step, as represented by the different box
sizes. c) After selection of a number of events using rejection-free kMC, the
final deformed state of the pillar reflects the operation of pre-existing
dislocations as well as a surface nucleation site. Arrow indicates most
probable event.

Then the next time step is computed as:

$\delta t^{n+1}=-\frac{log\xi_{2}~{}\Delta\varepsilon^{*}}{r_{tot}}$ (21)

such that the time is advanced as:

$t^{n+1}=t^{n}+\delta t^{n+1}$ (22)

$\xi_{1}$ and $\xi_{2}$ are random numbers uniformly distributed in $(0,1]$
and $\Delta\varepsilon^{*}$ is a normalization factor
($\Delta\varepsilon^{*}<=1$) to adjust the time step within physical time
scale bounds (Yu et al., 2021). Once an event is executed, the strain fields
and stress fields are updated.

##### Time scale bounds

The absolute time scale emanating from the sampling of eq. 17 using eqns. 18
through 21 represents the maximum time step compatible with the rate equation
under consideration. As shown by Yu et al. (2021), it is indeed possible to
use an arbitrarily smaller time step without invalidating the method. A
physically-reasonable upper bound for the time step is that of the total time
for a single dislocation to sweep the pillar section:

$t^{i}=\frac{D_{i}}{v^{i}}$ (23)

In addition, the existence of $t^{i}$ also poses a constraint. On each plastic
event, the eigenstrain applied cannot be larger than that of eq. 4, which
corresponds to a dislocation fully sweeping a pillar on a time $t^{i}$. Thus,
$t^{i}$ is an upper bound for the time increment $\delta t^{n+1}$. As
consequence, $\delta t^{n+1}/t^{i}~{}\leq 1$. In practice, the eigenstrain
applied once a plastic event is selected is:

$\boldsymbol{\varepsilon}^{EIG~{}(i)}(\mathbf{x})=\begin{cases}(\frac{\delta
t^{n+1}}{t^{i}})(\frac{b}{h})~{}(\hat{s}\otimes\hat{n}),~{}if~{}\mathbf{x}~{}\epsilon~{}plastic~{}region\\\
0,~{}elsewhere\end{cases}$ (24)

##### Weakest link principle and deformation localization

Under the weakest link principle, plastic flow occurs when the stress is high
enough to activate the weakest source available in the crystal (Norfleet et
al., 2008). In our model, weaker sources with a lower activation stress will
tend to produce higher dislocation velocities (eq. 13) and thus render higher
plastic displacement rates (eq. 12). In consequence, a sample with several
SASs, out of which one is significantly longer than the others, will result in
a higher probability of activation of such source, thus favoring localization
in the corresponding plastic region. In contrast, for a sample with several
SASs of commensurate length, the kMC selection process will translate into the
sequential activation of several SASs, favoring delocalization of deformation
into several plastic regions. These aspects are further discussed in Sec. 4
after the results presented in Sec. 3.

###  2.4 Solving the mechanical problem with an FFT algorithm

In order to obtain the strain and stress distribution in the pillar, which
drives the activation of slip in the plastic regions, an elastic problem with
eigenstrains has to be solved at each time step. An FFT-based approach will be
used for this purpose.

The total strain is the sum of the eigenstrains (input) and the elastic strain
$\boldsymbol{\varepsilon^{e}}$, which appears to enforce total strain
compatibility and stress equilibrium,

$\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon^{e}}+\boldsymbol{\varepsilon}^{EIG}.$
(25)

Under linear elasticity, stress is related to strain through

$\boldsymbol{\sigma}=\mathbb{C}:\boldsymbol{\varepsilon^{e}}=\mathbb{C}:(\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{EIG})$
(26)

where $\mathbb{C}$ is the stiffness tensor of the crystal.

The mechanical equilibrium corresponds to

$\nabla\cdot[\mathbb{C}:(\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{EIG})]=\boldsymbol{0}$
(27)

and rearranging terms, the result is a linear differential equation

$\nabla\cdot[\mathbb{C}:\boldsymbol{\varepsilon}]=\\\
-\nabla\cdot[\mathbb{C}:\boldsymbol{\varepsilon}^{EIG}].$ (28)

If the right-hand side of Eq. (28) is known, as in our case, the total strain
$\boldsymbol{\varepsilon}$ can be directly obtained by the convolution of this
term with the corresponding Green’s function derivative. Moreover, if the
medium is homogeneous and periodic, closed expressions for Green’s function
exist in Fourier space, and the solution is just a multiplication.

However, when considering the present problem in a FFT framework, the actual
geometry has to be embedded in a cuboidal periodic domain $\Omega$ which also
contains the outer free space. In this case, the domain can be considered as
an heterogeneous microstructure formed by two phases and this microstructure
can be described by

$\mathbb{C}=\begin{cases}\mathbb{C}_{1},~{}if~{}\mathbf{x}~{}in~{}\Omega_{1}\\\
\mathbb{C}_{2},~{}if~{}\mathbf{x}~{}in~{}\Omega_{2}\end{cases}$ (29)

In the simulation domain, the metallic pillar occupies the region
$(\Omega_{1})$ and is surrounded by an infinite compliant medium which
occupies the region $(\Omega_{2})$. For the solution of Eq. 28, considering
Eq. 29, we follow the Fourier-Galerkin approach. The details of the method can
be found in (Vondřejc et al., 2014; Zeman et al., 2017) and here only the
final equations are recalled. The problem is discretized in a box of
$[N_{x},N_{y},N_{z}]$ voxels which correspond to the same number of
frequencies in the discrete Fourier space. The value of $\mathbb{C}$, strain
and stress fields are approximated by trigonometrical polynomial, whose
coefficients are given by the FFT of the value of these fields at the voxels.
Let $\mathcal{G}$ be the projection operator that, by convolution, extracts
the compatible part of a tensor field (Vondřejc et al., 2014). The mechanical
equilibrium of the discrete stress field translates into,

$\mathcal{G}\ast\boldsymbol{\sigma}=\boldsymbol{0}$ (30)

where $\ast$ indicates a convolution. Introducing the value of
$\boldsymbol{\sigma}$ from eq. 26,

$\mathcal{G}\ast(\mathbb{C}(\mathbf{x}):(\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}^{EIG}))=\boldsymbol{0}.$
(31)

The total strain $\boldsymbol{\varepsilon}$ is decomposed into a homogeneous
average strain tensor, represented by the macro-scale applied strain
$\boldsymbol{E}^{T}$, and a periodic fluctuating micro-scale strain field
$\boldsymbol{\tilde{\varepsilon}}$, which is the unknown.

$\boldsymbol{\varepsilon}=\boldsymbol{\tilde{\varepsilon}}+\boldsymbol{E}^{T}$
(32)

Combining the last two equations results in a linear system in which the
strain fluctuation at the voxels are the unknowns.

$\mathcal{G}\ast(\mathbb{C}(\mathbf{x}):\boldsymbol{\tilde{\varepsilon}})=-\mathcal{G}\ast(\mathbb{C}(\mathbf{x}):(\boldsymbol{E}^{T}-\boldsymbol{\varepsilon}^{EIG}))$
(33)

This equation can be solved easily by transformation to Fourier space where
the operator has a closed expression and convolutions are transformed into
products. If a linear discrete operator $\mathcal{G}_{C}(\bullet)$ is defined
as

$\mathcal{G}_{C}(\bullet)=\mathcal{F}^{-1}(\hat{\mathcal{G}}:\mathcal{F}(\mathbb{C}(\mathbf{x}):(\bullet)))$
(34)

with $\mathcal{F}$ and $\mathcal{F}^{-1}$ the discrete Fourier transform
operations and $\hat{\mathcal{G}}$ the projection operator in Fourier space,
the equation to solve is

$\mathcal{G}_{C}(\tilde{\varepsilon})=-\mathcal{G}_{C}(\boldsymbol{E}^{T}-\boldsymbol{\varepsilon}^{EIG})).$
(35)

To solve this linear problem efficiently, an iterative solver has to be used.
The convergence rate and quality of the solution depend on the phase property
contrast, which here is infinite because one of the phases is empty. Following
Lucarini et al. (2022) this problem can be solved efficiently and with minimal
noise by using a discrete projection operator (here rotated discrete
$\mathcal{G}$-operator (Willot, 2015)) and a minimum residual iteration method
as linear solver.

It is important to highlight that the framework presented is thought to study
general geometrical domains in 3-D with heterogeneous microscopic fields.
While this may not be regarded as fundamental for a simple cylindrical pillar,
it is essential for the modeling of complex geometries, such as nanoporous
metals and nanoarchitected meta-materials.

For a set of initialization parameters, our method can be summarized by the
algorithm 1 presented in the Appendix.

##  3 Application to single-crystalline W pillar compression

Tungsten is a prototypical refractory metal with a high melting point and good
corrosion resistance. Its outstanding mechanical properties at high
temperatures, together with its low sputtering rate, make this material an
ideal candidate for plasma-facing components in fusion energy devices (Rieth
et al., 2013). In addition, tungsten is one of the most studied bcc metals by
micropillar compression testing (Schneider et al., 2009; Kim et al., 2010;
Abad et al., 2016; Srivastava et al., 2021) and its dislocation mobility laws
are well-documented (Cereceda et al., 2013, 2016; Po et al., 2016). These
aspects make W an ideal candidate to test the applicability of our framework.

Single crystalline tungsten pillars were simulated using the framework
described in the previous section. The pillar dimensions were varied in the
range from 25 nm to 1 $\mu$m at a strain-rate of $10^{-3}$ $s^{-1}$. The
crystalline orientation and loading axis were consistent with the [100]
direction. In all cases, an initial dislocation density of $\rho=5\cdot
10^{12}~{}m^{-2}$ was taken. Note that such an election renders an initial
dislocation count $N_{d}$ in the range of 1 for the smallest pillars to 8 for
the largest ones. The number of potential nucleation sites was taken as
$N_{SN}=128$. The parameters of the surface nucleation model are presented in
Table 1. The parameters of the physically-based dislocation mobility law are
presented in Table 2. A list of the slip systems considered can be found in
Table 3. The models used were discretized in a grid of $128\cdot 128\cdot 128$
voxels with a kMC normalization factor of $\Delta\varepsilon^{*}=10^{-3}$. See
Supplementary Material file for sensitivity of the results to the
discretization.

In particular, we focus on three important aspects: size effects (Sec.3.1),
strain-rate effects (Sec.3.2) and mechanical annealing effects (Sec.3.3). For
validation purposes, the results of the model are later compared with
experimental results in Sec. 4.

###  3.1 Size effects

Figure 3 presents the stress-strain curves for all the simulations performed
in this study. For each of the pillar dimensions probed, ten simulations were
conducted. Blue lines correspond to source truncation (SAS) dominated cases,
whereas red lines correspond to surface nucleation-dominated ones. The
envelope of the mechanical response for all the stress-strain curves is
illustrated by the corresponding colored areas. Black lines correspond to the
average of the results.

For the largest pillars, with a diameter of 1 $\mu m$, the plastic regime is
characterized by a rather smooth output with minimum scatter among the curves.
This is a consequence of having several single-arm sources distributed along
the sample, allowing the activation of several slip planes in a sequential
way. For 500 nm diameter pillars, the plastic regime is again smooth, but now
the curves appear more scattered, as judged by the upper and lower limits of
the blue-shaded region. For a smaller pillar diameter, the average number of
dislocations is close to one. As their pinning position is randomly selected,
significant scatter is then translated to the loading curves. The mean flow
stress is higher than before, as expected for a decreased pillar diameter.

For 200 nm diameter pillars, again a significant scatter is seen for the
loading curves. Interestingly, one case rendered a significantly high flow
stress, associated to surface nucleation of dislocations instead of the
activation of the SAS mechanism. This corresponds to the red curve, for which
a SAS length of 8 nm was randomly generated. The short length thus renders an
activation stress that is higher than that of the SN model. This behavior is
further analyzed in the Supplementary Material.

For 100 nm diameter pillars, now several cases appear both for the activation
of SA sources (blue) as well as surface nucleation (red). Both the mean flow
stresses for the SAS cases as well as the average (black) increase with
respect to the previous pillar diameter. For 50 nm diameter pillars, the
surface nucleation mechanism is dominant, with few cases of activation of SA
sources at higher stresses compared to previous dimensions. Finally, for 25 nm
diameter pillars, surface nucleation was found to be the only operating
mechanism. Note that as surface nucleation becomes the dominant mechanism, the
curves show smaller scatter. This is consistent with experimental observations
(Kiener and Minor, 2011; Huang et al., 2011).

Figure 3: Stress-strain response for W micro/nanopillars under compression.
For each of the pillar dimensions probed, ten simulations were conducted.
Shaded regions correspond to the envelope of the mechanical response for all
the stress-strain curves, where blue corresponds to single-arm source
dominated cases and red corresponds to surface nucleation dominated cases.
Blue and red lines correspond to the average of each subset (SAS-dominated and
SN-dominated, respectively), whereas black lines correspond to the average of
all the results, irrespectively of their nature.

###  3.2 Strain-rate effects

In order to evaluate the strain rate sensitivity of the flow strength, W
pillar simulations were conducted at different strain rates
($10^{-3}/s,~{}10^{-2}/s$ and $10^{-1}/s$) and for different pillar diameters
(1 $\mu$m, 500 nm, 200 nm and 100 nm). Figure 4 presents the logarithm of the
flow stress as function of the logarithmic applied strain rate obtained by the
model for pillars of different diameter and compared to that of Srivastava et
al. (2021) at room temperature. The log-log flow stress data was then fitted
with a power-law function of the strain rate,
$\sigma\propto\dot{\epsilon}^{m}$. The strain rate sensitivity parameter $m$
is then $m=\partial(ln\sigma)/\partial(ln\dot{\epsilon})$.

The resulting curves are approximately linear in the double-logarithm diagram,
indicating a strain rate sensitivity $m$ approximately constant in the strain
range explored. Comparison with Srivastava et al. (2021) shows not only
similar strain rate sensitivity parameters but also a similar transition trend
from positive to null values as the pillar diameter decreases.

Figure 4: Size-dependent strain-rate effects on the flow stress of W pillars.
Empty circles and dashed lines correspond to experimental results by
Srivastava et al. (2021). For comparison purposes, the colors are the same as
the ones used in Srivastava et al. (2021) and correspond to: blue - 1 $\mu$m;
red - 500 nm; magenta - 200 nm; green - 100 nm.

###  3.3 Mechanical annealing

In addition to single-arm sources, pre-existing pure screw dislocations can
also be considered, as presented in the methods section. These are valuable
for capturing mechanical annealing effects (Shan et al., 2008; Huang et al.,
2011; Shan, 2012) associated with dislocation starvation (Greer and Nix,
2006). Figure 5 presents three distinctive stress-strain curves corresponding
to 100 nm diameter pillar simulations where one pre-existing screw dislocation
has been considered. The red curve corresponds to a pillar without SASs and a
screw segment located in the middle of the cross-section (SN-dominated - case
1), while the green curve corresponds to a pillar without SASs and a screw
segment located halfway between the center of the pillar and its surface (SN-
dominated - case 2). In both cases, after the operation of the screw
dislocation, the source has exhausted. In consequence, the simulated test
continues elastically until activation of the SN mode. A similar scenario
would have been found for a case in which a pre-existing screw dislocation
habits a pillar with a SAS of sufficiently small length. In our
implementation, the position of the screw segment is randomly chosen in the
pillar cross section considering a uniform distribution. Specific knowledge of
the material under consideration would allow to modify such distribution
accordingly and without invalidating our methodology. Differences in load
drops for the activation of screw segments are discussed in Sec. 4. In
contrast, the blue curve corresponds to a pillar with a screw dislocation in
the middle of the cross section and a sufficiently long SA source. Upon
loading, and after an initial elastic stage, the screw dislocation is
activated -recall that the critical resolved shear stress is lower than that
of SASs -. After exhaustion of the screw source, the mechanical response
proceeds elastically until reaching a stress corresponding to the CRSS of the
SA source, which gets activated several times for the rest of the simulation.

Figure 5: Stress-strain response for 100 nm diameter pillar simulations
containing a pure screw dislocation. Red and green curves correspond to
pillars without SASs and a screw dislocation located in the middle of the
cross section (red - case 1) and halfway between the center of the pillar and
its surface (green - case 2). Blue curve corresponds to a pillar with a
sufficiently long SAS and a screw dislocation located in the middle of the
cross section.

##  4 Discussion

##### Comparison with other methods

Several aspects of our proposed framework are worth discussing. The framework
proposed aims at modelling the full elastic-plastic response with 3D
resolution of microfields. Its particular application to nano-pillars and
submicron-pillars allows to obtain the size effect in the yield like other
simple models (Parthasarathy et al., 2007; Jennings et al., 2013), but here
the microfields and post yield behavior are also recovered, with the aided
benefit of obtaining the 3-D deformation pattern (see below). Our model can
therefore be placed in between DDD approaches and continuum descriptions, such
as crystal plasticity. Regarding CP, continuum microfields are used but
deformation is localized here. In this aspect, our method shares a similarity
with the concept of thin discrete slip bands used by Wijnen et al. (2021) but
the nature of our model is stochastic and based on the displacement
discontinuity created by slip events. Like in DDD approaches (El-Awady et al.,
2009; Srivastava et al., 2013; Cui et al., 2014; Ryu et al., 2013, 2015,
2020), important aspects of the dislocation mobility laws and underlying
physics are considered. The selection of the plastic events does not take
place using the weakest-link principle, like in (El-Awady et al., 2009; Wijnen
et al., 2021), but a Monte Carlo selection process over the displacement rates
calculated using an approach similar to Orowan’s equation.

However, unlike DDD approaches, dislocations are not modeled explicitly, but
their effect is introduced by means of the concept of Eshelby inclusions, such
that the effect of dislocation gliding is introduced through an equivalent
eigenstrain. The total strain is thus the superposition of two fields,
($\boldsymbol{\varepsilon}$ and the inelastic strain obtained as a sum of
eigenstrains due to slip events $\boldsymbol{\varepsilon}^{EIG}$). Therefore,
the calculation of the strain and stress fields becomes linear, in contrast to
models based on localized crystal plasticity as (Wijnen et al., 2021). This
aspect, together with an inherently fast FFT solver, results in an efficient
and computationally cheap stochastic framework capable of simulating
experiments with 3D spatial resolution at strain rates lower than $10^{-3}/s$
on nano and submicron samples within a few hours on an off-the-shelf
workstation. The present implementation is based on small strains since we are
focusing on the initial flow of the pillars, but can be extended to finite
strains by assuming a non-linear elasticity energy and solving the FFT problem
in finite strains (Zeman et al., 2017). Other aspects for further extension of
the model include non-Schmid effects. This could be done following steps
already taken with recent advances on the development of a physically-informed
continuum crystal plasticity model for tantalum (Lee et al., 2023). The
framework would also allow for the incorporation of twinning as a
complementary deformation mechanism. Discrete dislocation plasticity models
have proven useful for the understanding of dislocation slip-mediated twinning
mechanisms (Wang et al., 2022) and could be used as input for the extension of
the framework presented here.

##### Flow stress size dependence

The proposed framework allows to obtain the statistical distribution of the
mechanical response of a pillar as a function of its size, as shown in Fig 3.
To analyze these results, in Figure 6 the flow stress for different pillar
sizes obtained with our model is compared to selected experimental references
(Schneider et al., 2009; Abad et al., 2016; Srivastava et al., 2021; Kim et
al., 2010) in a double logarithm representation. Flow stress values presented
correspond to the average stress in the plastic regime of the curves presented
in Fig. 3, while vertical bars correspond to standard deviation of all the
results presented in Fig. 3, capturing the dispersion of results produced by
the stochastic selection of the SAS length. Schneider et al. (2009) performed
load-controlled compression tests on [100]-oriented W pillars spanning from
the nm scale (200 nm) to the micrometer scale (6 $\mu m$). Later on, the same
group produced an experimental study probing temperature effects on W pillars
with similar orientation (Abad et al., 2016). Our model successfully captures
these experimental measurements as well as an overall $\sigma\propto
1/\lambda$ dependence of the yield strength in the 200 nm - 1 $\mu m$ range
(black dashed line). Srivastava et al. (2021) performed displacement-
controlled compression tests on [100]-oriented W nanopillars in the range of
100 nm - 1 $\mu m$. Their results are systematically lower than those reported
in Schneider et al. (2009); Abad et al. (2016) yet following a similar yield
strength trend increase with decreasing size ($\sigma\propto 1/\lambda$). Kim
et al. (2010) performed displacement-controlled compression tests on
[100]-oriented W nanopillars in the range of 200 - 900 nm. Their results
represent a lower bound for the experimental data set chosen. Yet their values
again follow similar trends ($\sigma\propto 1/\lambda$). To the best of the
authors knowledge, there are no reports on [100]-oriented W pillars with
diameters below 100 nm, probably due to the associated experimental
difficulties. At the limit of a 100 nm, we predict higher yield stresses than
those reported in Srivastava et al. (2021), yet our estimations agree with
their measurements if we take into account the standard deviation.

For diameters below 200 nm, the exponential dependence $\sigma\propto
1/\lambda$ yields to a surface nucleation-dominated trend (red dashed line),
indicating a crossover between predictions from the SAS and SN models. The
critical sample size observed is approximately 120 nm. Below this threshold,
the flow stress determined by the SN model is lower than that by the SAS
model, highlighting the comparative difficulty of activating a single-arm
source over nucleating a dislocation from a free surface, emphasizing a
preference for surface nucleation. In contrast, when dealing with samples that
exceed this critical size, plasticity is typically ruled by the activation of
pre-existing dislocation sources, such as those represented by the SAS
mechanism. A handful of previous studies had focused on the critical size for
bcc and fcc metals using both experimental and computational techniques. For
fcc materials, Hu et al. (2019) used DDD simulations to predict the the flow
stress and dominant yielding mechanisms in Cu nanopillars
($100~{}nm<D<800~{}nm$), reporting on a cross-over from SAS activation to SN
operation as D reaches $\approx 110nm$. Similarly, Shan et al. (2008) utilized
transmission electron microscopy (TEM) to investigate the behavior of pre-
existing dislocations in Ni pillars of different diameters. The authors
observed that in the case of a 160 nm diameter pillar, the pre-existing
dislocations gradually exited the pillar, followed by the emergence of new
dislocations, which can be attributed to the phenomenon known as SN.
Conversely, in the case of pillars with a diameter of 290 nm, dislocations
persisted even after undergoing compression, which can be tentatively
associated with SAS. For bcc metals, experiments conducted on bcc Mo pillars
suggest a critical size of around 200 nm (Huang et al., 2011; Shan, 2012). We
conclude that our critical size prediction for [100]-oriented W pillars is in
agreement with previous observations for other metals at these scales.

Figure 6: Log-log size-dependence of the flow stress compared with
experimental literature. Black dashed line indicates $1/\lambda$ dependence,
in agreement with a SAS dominated regime. Blue dashed line indicates the
estimated power law exponent, while light-blue shadow indicates $\pm$ values
($\alpha_{W}\approx 0.31\pm 0.11$). Red dashed line indicates a surface
nucleation (SN) trend. Its intersection indicates a cross-over transition
size. Vertical bars correspond to standard deviation.

The results can be interpreted through the well-known power law relation,
namely

$\sigma_{y}=\sigma_{0}+k_{1}D^{-\alpha_{W}}$ (36)

where $\sigma_{0}$ is the strength of the bulk material, $k_{1}$ is a constant
and $\alpha_{W}$ is the power law exponent. In practice, $\sigma_{0}$ is
usually neglected (Brinckmann et al., 2008; Abad et al., 2016; Srivastava et
al., 2021) as the second term in eq. 36 dominates at these small scales. As a
reference, Srivastava et al. (2021) report on exponent of 0.319 for their
$10^{-3}~{}s^{-1}$ experiments, close to the $\alpha_{W}=0.33$ value obtained
by Abad et al. (2016) for the same strain rate. In turn, Kim et al. (2010)
informed a power law exponent in the range of $0.33<\alpha_{W}<0.55$. On a
lower end, Schneider et al. (2009) report on $\alpha_{1}=0.21$. By linear
fitting eq. 36 on the SAS-dominated regime ($D>120~{}nm$), we obtain an
exponent $\alpha_{W}\approx 0.31\pm 0.11$, in agreement with most of the
experimental evidence. This further supports the ability of the single-arm
source model to understand the power-law relation on the strength vs size
dependence of micropillars (Lee and Nix, 2012). The red dashed line,
connecting the average yield stress values obtained for the smaller pillars,
presents a slope that seems to be an indication of a possibly weaker size
dependence.

The underlying cause for the values of $\alpha_{W}$ are worth a discussion.
Brinckmann et al. (2008) used the screw dislocation cross-slip mechanism
concept, derived from high strain-rate MD and DD simulations (Weinberger and
Cai, 2008; Greer et al., 2008), to provide an explanation for their bcc Mo
pillar compression test results. They argued that the relatively low power law
exponent of bcc metals ($\alpha_{Mo}=0.45$) could be attributed to such
effects. Interestingly, Kim et al. (2010) performed pillar compression tests
on W, reporting on strikingly similar power-law exponents ($\alpha_{W}=0.44\pm
0.11$). One could then assume that the same mechanism takes place in bcc W.
Yet, our power-law exponent estimations are on the lower bound of this regime,
while our model does not incorporate a single characteristic of the cross-slip
mechanism presented in (Weinberger and Cai, 2008). Therefore, our results
suggest that screw dislocation mobility alone could indeed provide adequate
power-law exponents provided the underlying physics of the dislocation
mobility law are fully incorporated ( eq. 13).

It is worth highlighting that previous experimental observations on bcc Mo
suggest that size effect itself has a strong size effect (Huang et al., 2011;
Shan, 2012). In their in-situ TEM bcc Mo studies, Huang et al. (2011) found a
strong increase in the power law exponent as the pillar diameter decreases
beyond 200 nm. In contrast, we predict a break down in the power-law exponent
as the W pillar diameter decreases below 200 nm. This difference in behavior
could be attributed to different activation energies between Mo and W. Small
variations in activation energy can strongly influence yield stress at the
nanoscale. A break down in exponent can also be inferred from dislocation
dynamics simulations on Cu (Hu et al., 2019). We attribute the break down in
power law exponent to the change of mechanism, from SAS-dominated to surface
nucleation dominated, for which a weaker stress dependence on diameter is
expected (Zhu et al., 2008; Huang et al., 2015; Hu et al., 2019). Certainly,
this alternative explanation of the results leaves space for future research
and informed discussion.

##### Strain-rate sensitivity

Our results show a size dependence of the strain rate sensitivity parameter,
which changes from positive $m$ to nearly null values (Figure 4). Huang et al.
(2015) and Srivastava et al. (2021) observed the same positive to null trends
for their bcc Fe and W pillar experiments, respectively. The strain rate
sensitivity exponent $m$ depends on the activation volume $V_{a}$ and on the
resolved stress on the active plane $\tau$ through (Kocks et al., 1975):

$m=\frac{k_{B}T}{\tau V_{a}}.$ (37)

In the previous equation, both the critical shear $\tau$ and the activation
volume $V_{a}$ can be size dependent, and the origin of size dependency of $m$
might have both contributions. Regarding the effect of activation volume,
Huang et al. (2015) rationalized that in Fe this might be the main
contribution to size effect in $m$. This size dependent $m$ in Fe arises from
a competition between a size-independent friction term -that rules for large
pillar diameters (D $>$ 500 nm)- and a size-dependent term -that rules for
smaller pillar diameters-, with larger activation volumes for increasingly
smaller pillars (100 nm $<$ D $<$ 500 nm). They also argue that as the pillar
diameter further decreases (D $<$ 100 nm), the activation volume becomes again
smaller. Still, Fe is known to have an activation volume that may display
larger variations with stress compared to other bcc metals with higher melting
point (Christian and Masters, 1964). Tungsten, in constrast, and particularly
for low homologous temperatures, displays a very low activation volume
($V_{a}\leq 10b^{3}$) (Kiener et al., 2019), consistent with the dominance of
the kink mechanism at low temperatures. Our calculations of the activation
volume (see Suppl. Mat.) render that for the kink-pair activated mechanism,
the stress-dependent activation volume is in the range of 4 $b^{3}$ $<$
$V_{a}^{kink}$ $<$ 11 $b^{3}$ for the stress range of the SAS-dominated
regime, whereas for the surface nucleation mechanism is in the range of 2.5
$b^{3}$ $<$ $V_{a}^{SN}$ $<$ 4.5 $b^{3}$ for the stress range of the SN-
dominated regime.

Therefore, considering the relatively low variations in the activation volume
$V_{a}$ in W, changes in the strain rate sensitivity exponent $m$ can be
understood based on the size-dependence of the stress reached on the pillar
due to the size dependent strength of SAS, as well as on the stress-dependence
of the dislocation mobility law. Similar conclusions on the importance of
$\tau$ were reached by Srivastava et al. (Srivastava et al., 2021) for their
strain-rate sensitivity studies on W pillars.

Future investigations using our proposed framework could also focus on
temperature effects, as a recently developed CP model using similar
dislocation mobility laws has proven useful in exploring the temperature
dependence of deformation localization in irradiated tungsten (Li et al.,
2021).

##### Mechanical annealing

Earlier MD simulation studies predict that such behavior cannot take place in
bcc pillars. Such simulations suggested that the combined effects of image
forces, together with the dislocation core structure, would facilitate
dislocation multiplication by a surface-assisted cross-slip mechanism
(Weinberger and Cai, 2008). SEM observations of multiple slip traces on
compressed bcc pillars seem to support this mechanism (Schneider et al.,
2009). On the other hand, in-situ TEM compression of Mo nanopillars (Huang et
al., 2011; Shan, 2012) and Fe–3% Si pillars (Zhang et al., 2012) show clear
evidence that mechanical annealing by dislocation starvation can indeed take
place in bcc metals. Our proposed computational framework takes into account
such scenarios. Interestingly, load drops due to mechanical annealing are, in
general, less marked than in our Figure 5, as shown in (Shan et al., 2008).
Other works do report stress drops of the order of 0.5 - 1.0 GPa (Huang et
al., 2011; Zhang et al., 2012). Both situations can be taken into account in
our method by an adequate selection of the statistical distribution of pre-
existing screw segments.

##### Deformation patterns and localization

Figure 7 presents the typical deformation patterns obtained in our
simulations, together with a comparison with a selected reference (Srivastava
et al., 2021). For our largest pillar diameter ($D=1\mu m$), the deformation
pattern reveals the activation of several SASs on a variety of slip systems
(see Supplementary Material movie), in agreement with scanning electron
microscopy (SEM) analysis of W micron-sized pillars (Schneider et al., 2009;
Abad et al., 2016). As the pillar diameter decreases to the submicron regime
($D=500~{}nm$ and $D=200~{}nm$) and for the dislocation density chosen, only
one SAS is present. Upon activation, significant deformation takes place and
concentrates on a single slip band, in agreement with SEM micrographs of W
pillars with comparable diameters (Srivastava et al., 2021). Further decrease
of pillar diameter into the nanometer regime ($D=100~{}nm$) triggers surface
nucleation instead of SAS activation. The former takes place on multiple slip
systems.

As explained in Sec. 2.3, our model favors the activation of weaker sources.
Localization of plastic activity in one or a few number of planes depends on
the number of potential slip sites (e.g. number of SAS) and the statistical
distribution of the slip resistance on that planes (linked with the length of
each SAS). Our 1 $\mu m$ diameter pillars with $\rho=5\cdot 10^{12}~{}m^{-2}$
have multiple SASs of variable length. They show a tendency for a rather
delocalized (homogeneous) distribution of deformation. In contrast, 200 nm -
500 nm diameter pillars with similar dislocation density have 1 - 2 SASs on
average. In such cases, localization is favored (Figure 7). Cui et al. (2018)
presented a 2D Monte-Carlo model of DD source activation coupled with
crossslip channel widening to reproduce and physically explain the transition
in the mechanism of plastic flow localization in irradiated materials from
irradiation-controlled to dislocation source-controlled. They showed that as
the size decreases, the spatial correlation of plastic deformation decreases
due to weaker dislocation interactions and less frequent cross-slip, thus
producing thinner dislocation channels. For a low irradiation damage, they
show that as the diameter decreases from 1.5 $\mu m$ to 300 nm, deformation
transitions from rather homogeneous to highly localized. Our results broadly
agree with this picture.

##### Statistical analysis of burst displacement

In order to analyze the statistics of the plastic bursts obtained in the
simulations, the complementary cumulative distribution function (CCDF) of the
burst displacement magnitude $\Delta U$ is obtained. The CDDF is computed
using the method presented by Cui et al. (2016) and its result for a 1 $\mu$m
diameter pillar is represented in Figure 8. The Figure suggest that, under the
displacement-controlled conditions (strain control) explored here, the system
does not exhibit power-law scaling, in agreement with Cui et al. (2016). In
addition, the data spans less than two orders of magnitude. In contrast,
typical load-displacement controlled compression experiments usually span
several orders of magnitude, with bursts often exceeding 10 nm (Alcalá et al.,
2020; Srivastava et al., 2021). Detailed tracking of the plastic strain rate
(Inset of Figure 8) shows that the dynamical behavior is quasiperiodic, in
agreement with displacement-controlled simulations (Cui et al., 2016) and
experiments (Papanikolaou et al., 2012). In contrast, experiments under load-
controlled conditions lead to a completely different scenario, where
plasticity does not lead to marked stress drops, allowing to trigger the
operation of several SAS simultaneously and thus exhibiting power law scaling.
As a result, avalanche events can develop, leading to a highly correlated
dynamical response, consistent with the concept of self-organized criticality
(Bak et al., 1988). The framework presented here, inherently implies the
activation of one source at a time, leading to a quasiperiodic response that
it is consistent with displacement-control conditions.

##### On the incorporation of strain hardening effects

In all the cases presented in Figure 3, the flow stress keeps constant, in
agreement with the assumptions introduced in Sec. 2. This is consistent with
an scenario in which the slip planes swept by the SAS exhibit a scarcity of
dislocations in intersecting planes, and where the pinning point of the SAS is
stable under the prevailing stress levels. Consequently, the activated SAS
does not encounter other defects, hence preventing the formation of new
junctions or the reduction of the source length. In addition, experiments on W
single crystalline pillars (Schneider et al., 2009) do not show hardening for
the relatively low strains explored in Figure 3.

For larger samples, our method could be extended to incorporate hardening
effects by including a Taylor hardening-like term in eq. 6 and a dislocation
density evolution law (Cui et al., 2014). The latter should incorporate size-
dependent terms and bulk-like terms. Size-dependent terms include the
generation of dislocations due to the operation of existing sources ($\propto
1/\lambda$) and the escape of dislocations reaching a free surface
($\propto-1/D$) (Cui et al., 2014). Bulk-like terms include dislocation
multiplication due to forest dislocations ($\propto\sqrt{\rho}$) and the
annihilation of closely spaced dislocations of opposite signs ($\propto-\rho$)
(Devincre et al., 2008). If self and pair interactions must be taken into
account (e.g. for micron scale pillars above 1 $\mu$m), the corresponding
hardening coefficient matrix can also be incorporated (Yu et al., 2021; Lee et
al., 2023). It must be emphasized that for submicron pillars, simple
calculations show that the contributions of the bulk-like terms are nearly an
order of magnitude smaller than the size-dependent terms (Cui et al., 2014).
Thus, by neglecting the former and provided the dislocation generation rate is
balanced with the escape rate from free surface which, together with a
relatively low dislocation density, the results presented in Sec.3.1 would be
reproduced.

To wrap up, the results presented here focus solely on tungsten and a
relatively simple geometry. Yet, these are not restrictions. For instance,
FeCrAl alloys and refractory multi-principal element alloys (RMPEAs) are also
promising candidate materials for use under extreme conditions in the nuclear
industry. FeCrAl alloys have already been probed using discrete dislocation
dynamics (Pachaury et al., 2023) as well as with crystal plasticity approaches
using Arrhenius-type rate equations (Gong et al., 2023). With respect to
RMPEAs, kink-migration of screw dislocations have been shown to play a
dominant role in their deformation (Balbus et al., 2024) and it is expected
that developments in dislocation mobility laws (Shen and Spearot, 2021) would
allow for systematic investigations on these promising alloys using upper
scale methods. Our proposed framework offers flexibility to explore the
micropillar compression behavior of such alloys by the incorporation of
specific dislocation mobility laws. To further extend its potential, the use
of a 3D FFT solver offers opportunities for the exploration of a variety of
complex geometries, including nanoporous metals and nanoarchitected meta-
materials.

Figure 7: Residual deformation patterns for different pillar diameters
explored, together with SEM observations by Srivastava et al. (2021).
Reprinted with permission from Elsevier. Frame colors correspond to dominating
deformation mechanism (Blue - Single-Arm source, Red - Surface nucleation).
Figure 8: Complementary cumulative distribution function (CCDF) of burst
displacement for one example of our 1 $\mu$m diameter pillars. Inset,
corresponding evolution of plastic strain rate showing typical quasiperiodic
strain bursts, consistent with displacement-controlled compression conditions.
Note the breakdown from power law scaling (dashed line).

##  5 Conclusions

A novel approach is presented to simulate the deformation of submicron
specimens, accounting for the discrete events produced by the slip of internal
dislocation and surface nucleated ones in a stochastic manner. The framework
considers the slip events as eigenstrain fields that produce a displacement
jump across a slip plane and whose activation is driven by a Monte Carlo (MC)
method. Physically-based laws are incorporated to account for activation
probabilities, dislocation mobility and surface nucleation.

Two factors of stochasticity are taken into account: a random selection of the
position of defects and a random selection of plastic events after a MC
process on a sampling array of possible displacement rates, computed using a
similar concept as in Orowan’s equation. Implementation on a fast FFT solver,
results on an efficient, computationally-cheap algorithm which allows to
simulate accurately the deformation of a specimen in a wide range of strain
rates and sizes in a fraction of the time needed using techniques as MD or
DDD.

The framework developed has been applied to the compression of tungsten (W) on
[100]-oriented pillars in the 25 nm to 1 $\mu m$ diameter range. We find an
adequate agreement between simulations and experimental data. The study a
size-dependence of flow stress, characterized by distinct power-law exponents
for regimes dominated by source truncation (D $>$ 120 nm, $\alpha_{W}\approx
0.31\pm 0.11$) and a break down of the power law scaling as surface surface
nucleation becomes dominant (D $<$ 120 nm). Pre-existing dislocations
naturally compete with surface nucleation of new dislocations. The former are
favored for pillars above 120 nm diameter, whereas the latter is favored below
such critical size. In addition, strain-rate sensitivity effects are
adequately captured, including a size-dependence of the strain-rate
sensitivity exponent. This is attributed not only to changes in activation
volume and source strength but also to the stress-dependent dislocation
mobility law, whose form considers both the kink-pair regime, dominating at
the relatively low stresses found in submicron pillars, as well as the phonon
drag regime, that rules for very high stresses found in pillars of 200 nm and
below. The framework inherently implies the activation of one source at a
time, leading to a quasiperiodic response that it is consistent with
displacement-control conditions.

Finally, the framework presented here opens the possibility to provide
valuable analysis and interpretations of other small-scale testing techniques,
such as microtensile testing and microbending. Future developments of
dislocation mobility laws for FeCrAl alloys and multi-principal element alloys
would allow for a systematic investigation of these relevant metals, while the
possibility to include other micro and nanoscale geometries further extends
the potential applications map to include nanolattices and nanoporous
materials, among other nanostructures of interest.

##  6 Acknowledgements

This project has received funding from the European Union’s Horizon Europe
research and innovation programme under the Marie Sklodowska-Curie grant
agreement no. 101062254. Funded by the European Union. Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of the European Union. Neither the European Union nor the
granting authority can be held responsible for them.

##  Appendix A Appendix Section

Algorithm 1 kMC discrete slip approach with eigenstrains

1:procedure Initialization

  * 1.

Define pillar geometry, $\mathbb{C}(\mathbf{x})$, in the discrete domain.

  * 2.

Set positions, eigenstrain fields and strength of each plastic region $i$

  * 3.

Set loading conditions $\mathbf{E}^{T}(t)$,$\Delta\varepsilon^{*}$.

2:procedure Solution

3: while $t<t_{TOT}$ do:

4: Get rnd no. $\xi_{1}$, $\xi_{2}$ $\in$ (0,1]

5: Solve the equilibrium (Eq. 35) using FFT solver to obtain
$\boldsymbol{\sigma}$ (Eq. 26) , $|\tau^{i}|$ (Eq. 16)

6: for each region do

7: if $|\tau_{i}|>\tau_{CRSS}$ then

8: Compute: $P$ (SN sites only), $v^{i}$ (Eq. 13) , $t^{i}$ (Eq. 23) and
$\dot{u}^{i}_{p}$ (Eq. 12).

9: else

10: $v^{i}=0$

11: Compute accumulated rates $r_{j}$, $r_{tot}$

12: $\delta t^{n+1}=-\frac{log\xi_{2}~{}\Delta\varepsilon^{*}}{r_{tot}}$

13: if $\xi_{1}<\frac{r_{0}}{r_{t}}$ then

14: $\delta t^{n+1}=dt^{0}$

15: next event is elastic: $\boldsymbol{\varepsilon}^{EIG~{}(i)}=0$

16: else

17: next event is plastic

18: Compute eigenstrain: $\boldsymbol{\varepsilon}^{EIG~{}(i)}$ (Eq. 24)

19:

20: Update total eigenstrain field $\boldsymbol{\varepsilon}^{EIG}$ (Eq. 5)

The slip systems considered in this study are presented in Table 3.

$i$ | Slip system | $s^{i}$ | $n^{i}$
---|---|---|---
1 | $[1\bar{1}1](011)$ | $[1\bar{1}1]$ | $[011]$
2 | $[\bar{1}\bar{1}1](011)$ | $[\bar{1}\bar{1}1]$ | $[011]$
3 | $[111](0\bar{1}1)$ | $[111]$ | $[0\bar{1}1]$
4 | $[\bar{1}11](0\bar{1}1)$ | $[\bar{1}11]$ | $[0\bar{1}1]$
5 | $[\bar{1}11](101)$ | $[\bar{1}11]$ | $[101]$
6 | $[\bar{1}\bar{1}1](101)$ | $[\bar{1}\bar{1}1]$ | $[101]$
7 | $[111](\bar{1}01)$ | $[111]$ | $[\bar{1}01]$
8 | $[1\bar{1}1](\bar{1}01)$ | $[1\bar{1}1]$ | $[\bar{1}01]$
9 | $[\bar{1}11](110)$ | $[\bar{1}11]$ | $[110]$
10 | $[\bar{1}1\bar{1}](110)$ | $[\bar{1}1\bar{1}]$ | $[110]$
11 | $[111](\bar{1}10)$ | $[111]$ | $[\bar{1}10]$
12 | $[11\bar{1}](\bar{1}10)$ | $[11\bar{1}]$ | $[\bar{1}10]$
Table 3: Slip systems considered in this study

## References

  * Abad et al. (2016) Abad, O.T., Wheeler, J.M., Michler, J., Schneider, A.S., Arzt, E., 2016.  Temperature-dependent size effects on the strength of ta and w micropillars.  Acta Materialia 103, 483–494.
  * Alcalá et al. (2020) Alcalá, J., Očenášek, J., Varillas, J., A. El-Awady, J., Wheeler, J.M., Michler, J., 2020\.  Statistics of dislocation avalanches in fcc and bcc metals: dislocation mechanisms and mean swept distances across microsample sizes and temperatures.  Scientific reports 10, 19024\.
  * Bak et al. (1988) Bak, P., Tang, C., Wiesenfeld, K., 1988.  Self-organized criticality.  Physical review A 38, 364\.
  * Balbus et al. (2024) Balbus, G., Rao, S., Senkov, O., Payton, E., 2024\.  Orientation dependent plasticity of the refractory multi-principal element alloy monbti investigated via micropillar compression.  Acta Materialia 262, 119401\.
  * Brinckmann et al. (2008) Brinckmann, S., Kim, J.Y., Greer, J.R., 2008.  Fundamental differences in mechanical behavior between two types of crystals at the nanoscale.  Physical Review Letters 100, 155502\.
  * Capolungo and Taupin (2019) Capolungo, L., Taupin, V., 2019\.  Gd 3: generalized discrete defect dynamics.  Materials Theory 3, 1–21.
  * Cereceda et al. (2016) Cereceda, D., Diehl, M., Roters, F., Raabe, D., Perlado, J.M., Marian, J., 2016\.  Unraveling the temperature dependence of the yield strength in single-crystal tungsten using atomistically-informed crystal plasticity calculations.  International Journal of Plasticity 78, 242–265.
  * Cereceda et al. (2013) Cereceda, D., Stukowski, A., Gilbert, M., Queyreau, S., Ventelon, L., Marinica, M.C., Perlado, J., Marian, J., 2013\.  Assessment of interatomic potentials for atomistic analysis of static and dynamic properties of screw dislocations in w.  Journal of Physics: Condensed Matter 25, 085702.
  * Chen et al. (2015) Chen, L.Y., He, M.r., Shin, J., Richter, G., Gianola, D.S., 2015.  Measuring surface dislocation nucleation in defect-scarce nanostructures.  Nature materials 14, 707–713.
  * Christian (1983) Christian, J., 1983.  Some surprising features of the plastic deformation of body-centered cubic metals and alloys.  Metallurgical transactions A 14, 1237–1256.
  * Christian and Masters (1964) Christian, J.W., Masters, B., 1964\.  Low-temperature deformation of body-centred cubic metals. ii. mechanism of thermally activated flow.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 281, 240–257.
  * Córdoba et al. (2017) Córdoba, R., Lorenzoni, M., Pablo-Navarro, J., Magén, C., Pérez-Murano, F., De Teresa, J.M., 2017\.  Suspended tungsten-based nanowires with enhanced mechanical properties grown by focused ion beam induced deposition.  Nanotechnology 28, 445301\.
  * Cruzado et al. (2015) Cruzado, A., Gan, B., Jiménez, M., Barba, D., Ostolaza, K., Linaza, A., Molina-Aldareguia, J., Llorca, J., Segurado, J., 2015.  Multiscale modeling of the mechanical behavior of in718 superalloy based on micropillar compression and computational homogenization.  Acta Materialia 98, 242–253.
  * Cui et al. (2014) Cui, Y., Lin, P., Liu, Z., Zhuang, Z., 2014.  Theoretical and numerical investigations of single arm dislocation source controlled plastic flow in fcc micropillars.  International Journal of Plasticity 55, 279–292.
  * Cui et al. (2016) Cui, Y., Po, G., Ghoniem, N., 2016.  Controlling strain bursts and avalanches at the nano-to micrometer scale.  Physical review letters 117, 155502\.
  * Cui et al. (2018) Cui, Y., Po, G., Ghoniem, N., 2018.  Size-tuned plastic flow localization in irradiated materials at the submicron scale.  Physical review letters 120, 215501\.
  * Cui et al. (2015) Cui, Y.n., Liu, Z.l., Zhuang, Z., 2015.  Theoretical and numerical investigations on confined plasticity in micropillars.  Journal of the Mechanics and Physics of Solids 76, 127–143.
  * Dehm et al. (2018) Dehm, G., Jaya, B.N., Raghavan, R., Kirchlechner, C., 2018\.  Overview on micro-and nanomechanical testing: New insights in interface plasticity and fracture at small length scales.  Acta Materialia 142, 248–282.
  * Devincre et al. (2008) Devincre, B., Hoc, T., Kubin, L., 2008.  Dislocation mean free paths and strain hardening of crystals.  Science 320, 1745–1748.
  * El-Awady et al. (2009) El-Awady, J.A., Wen, M., Ghoniem, N.M., 2009.  The role of the weakest-link mechanism in controlling the plasticity of micropillars.  Journal of the Mechanics and Physics of Solids 57, 32–50.
  * Fan et al. (2021) Fan, H., Wang, Q., El-Awady, J.A., Raabe, D., Zaiser, M., 2021.  Strain rate dependency of dislocation plasticity.  Nature communications 12, 1845\.
  * Gong et al. (2023) Gong, M., Xie, D., Sun, T., Zhang, X., Shao, L., Wang, J., 2023.  Bridging microscale to macroscale mechanical property measurements of fecral alloys by crystal plasticity modeling.  International Journal of Plasticity 165, 103608.
  * Greer and De Hosson (2011) Greer, J.R., De Hosson, J.T.M., 2011\.  Plasticity in small-sized metallic systems: Intrinsic versus extrinsic size effect.  Progress in Materials Science 56, 654–724.
  * Greer and Nix (2006) Greer, J.R., Nix, W.D., 2006\.  Nanoscale gold pillars strengthened through dislocation starvation.  Physical Review B 73, 245410\.
  * Greer et al. (2008) Greer, J.R., Weinberger, C.R., Cai, W., 2008.  Comparing the strength of fcc and bcc sub-micrometer pillars: Compression experiments and dislocation dynamics simulations.  Materials Science and Engineering: A 493, 21–25.
  * Gu et al. (2021) Gu, Y., Eastman, D.W., Hemker, K.J., El-Awady, J.A., 2021\.  A statistical model for predicting size effects on the yield strength in dislocation-mediated crystal plasticity.  Journal of the Mechanics and Physics of Solids 147, 104245.
  * Healy and Ackland (2014) Healy, C.J., Ackland, G.J., 2014\.  Molecular dynamics simulations of compression–tension asymmetry in plasticity of fe nanopillars.  Acta Materialia 70, 105–112.
  * Hu et al. (2017) Hu, J., Liu, Z., Van der Giessen, E., Zhuang, Z., 2017.  Strain rate effects on the plastic flow in submicron copper pillars: Considering the influence of sample size and dislocation nucleation.  Extreme Mechanics Letters 17, 33–37.
  * Hu et al. (2019) Hu, J., Song, H., Liu, Z., Zhuang, Z., Liu, X., Sandfeld, S., 2019.  Predicting the flow stress and dominant yielding mechanisms: analytical models based on discrete dislocation plasticity.  Scientific Reports 9, 20422\.
  * Huang et al. (2011) Huang, L., Li, Q.J., Shan, Z.W., Li, J., Sun, J., Ma, E., 2011.  A new regime for mechanical annealing and strong sample-size strengthening in body centred cubic molybdenum.  Nature communications 2, 547\.
  * Huang et al. (2015) Huang, R., Li, Q.J., Wang, Z.J., Huang, L., Li, J., Ma, E., Shan, Z.W., 2015.  Flow stress in submicron bcc iron single crystals: sample-size-dependent strain-rate sensitivity and rate-dependent size strengthening.  Materials Research Letters 3, 121–127.
  * Jennings et al. (2013) Jennings, A.T., Weinberger, C.R., Lee, S.W., Aitken, Z.H., Meza, L., Greer, J.R., 2013\.  Modeling dislocation nucleation strengths in pristine metallic nanowires under experimental conditions.  Acta materialia 61, 2244–2259.
  * Jérusalem et al. (2012) Jérusalem, A., Fernández, A., Kunz, A., Greer, J.R., 2012\.  Continuum modeling of dislocation starvation and subsequent nucleation in nano-pillar compressions.  Scripta Materialia 66, 93–96.
  * Kiener et al. (2019) Kiener, D., Fritz, R., Alfreider, M., Leitner, A., Pippan, R., Maier-Kiener, V., 2019\.  Rate limiting deformation mechanisms of bcc metals in confined volumes.  Acta Materialia 166, 687–701.
  * Kiener and Minor (2011) Kiener, D., Minor, A.M., 2011\.  Source-controlled yield and hardening of cu (1 0 0) studied by in situ transmission electron microscopy.  Acta Materialia 59, 1328–1337.
  * Kiener et al. (2009) Kiener, D., Motz, C., Dehm, G., 2009.  Micro-compression testing: A critical discussion of experimental constraints.  Materials Science and Engineering: A 505, 79–87.
  * Kim et al. (2010) Kim, J.Y., Jang, D., Greer, J.R., 2010.  Tensile and compressive behavior of tungsten, molybdenum, tantalum and niobium at the nanoscale.  Acta Materialia 58, 2355–2363.
  * Kocks et al. (1975) Kocks, U.F., AS, A., MF, A., 1975.  Thermodynamics and kinetics of slip .
  * Konstantinidis et al. (2014) Konstantinidis, A.A., Aifantis, K.E., De Hosson, J.T.M., 2014.  Capturing the stochastic mechanical behavior of micro and nanopillars.  Materials Science and Engineering: A 597, 89–94.
  * Lee et al. (2023) Lee, S., Cho, H., Bronkhorst, C.A., Pokharel, R., Brown, D.W., Clausen, B., Vogel, S.C., Anghel, V., Gray III, G.T., Mayeur, J.R., 2023\.  Deformation, dislocation evolution and the non-schmid effect in body-centered-cubic single-and polycrystal tantalum.  International Journal of Plasticity 163, 103529.
  * Lee and Nix (2012) Lee, S.W., Nix, W.D., 2012.  Size dependence of the yield strength of fcc and bcc metallic micropillars with diameters of a few micrometers.  Philosophical Magazine 92, 1238–1260.
  * Li et al. (2021) Li, Z., Liu, Z., Zhuang, Z., Cui, Y., 2021.  Temperature dependent deformation localization in irradiated tungsten.  International Journal of Plasticity 146, 103077.
  * Lin et al. (2015) Lin, P., Liu, Z., Cui, Y., Zhuang, Z., 2015.  A stochastic crystal plasticity model with size-dependent and intermittent strain bursts characteristics at micron scale.  International Journal of Solids and Structures 69, 267–276.
  * Lin et al. (2016) Lin, P., Liu, Z., Zhuang, Z., 2016.  Numerical study of the size-dependent deformation morphology in micropillar compressions by a dislocation-based crystal plasticity model.  International Journal of Plasticity 87, 32–47.
  * Lucarini et al. (2022) Lucarini, S., Cobian, L., Voitus, A., Segurado, J., 2022\.  Adaptation and validation of fft methods for homogenization of lattice based materials.  Computer Methods in Applied Mechanics and Engineering 388, 114223.
  * Marichal et al. (2013) Marichal, C., Van Swygenhoven, H., Van Petegem, S., Borca, C., 2013\.  $\\{$110$\\}$ slip with $\\{$112$\\}$ slip traces in bcc tungsten.  Scientific reports 3, 2547\.
  * Mura (2013) Mura, T., 2013.  Micromechanics of defects in solids.  Springer Science & Business Media.
  * Ng and Ngan (2008) Ng, K., Ngan, A., 2008.  Stochastic nature of plasticity of aluminum micro-pillars.  Acta Materialia 56, 1712–1720.
  * Norfleet et al. (2008) Norfleet, D., Dimiduk, D., Polasik, S., Uchic, M., Mills, M., 2008.  Dislocation structures and their relationship to strength in deformed nickel microcrystals.  Acta Materialia 56, 2988–3001.
  * Pachaury et al. (2023) Pachaury, Y., Warren, G., Wharry, J.P., Po, G., El-Azab, A., 2023.  Plasticity in irradiated fecral nanopillars investigated using discrete dislocation dynamics.  International Journal of Plasticity , 103676\.
  * Papanikolaou et al. (2012) Papanikolaou, S., Dimiduk, D.M., Choi, W., Sethna, J.P., Uchic, M.D., Woodward, C.F., Zapperi, S., 2012.  Quasi-periodic events in crystal plasticity and the self-organized avalanche oscillator.  Nature 490, 517–521.
  * Parthasarathy et al. (2007) Parthasarathy, T.A., Rao, S.I., Dimiduk, D.M., Uchic, M.D., Trinkle, D.R., 2007.  Contribution to size effect of yield strength from the stochastics of dislocation source lengths in finite samples.  Scripta Materialia 56, 313–316.
  * Po et al. (2016) Po, G., Cui, Y., Rivera, D., Cereceda, D., Swinburne, T.D., Marian, J., Ghoniem, N., 2016\.  A phenomenological dislocation mobility law for bcc metals.  Acta Materialia 119, 123–135.
  * Raabe et al. (2007) Raabe, D., Ma, D., Roters, F., 2007.  Effects of initial orientation, sample geometry and friction on anisotropy and crystallographic orientation changes in single crystal microcompression deformation: A crystal plasticity finite element study.  Acta Materialia 55, 4567–4583.
  * Rao et al. (2008) Rao, S.I., Dimiduk, D., Parthasarathy, T.A., Uchic, M., Tang, M., Woodward, C., 2008\.  Athermal mechanisms of size-dependent crystal flow gleaned from three-dimensional discrete dislocation simulations.  Acta Materialia 56, 3245–3259.
  * Rao et al. (2007) Rao, S.I., Dimiduk, D., Tang, M., Uchic, M., Parthasarathy, T., Woodward, C., 2007\.  Estimating the strength of single-ended dislocation sources in micron-sized single crystals.  Philosophical Magazine 87, 4777–4794.
  * Rieth et al. (2013) Rieth, M., Dudarev, S.L., De Vicente, S.G., Aktaa, J., Ahlgren, T., Antusch, S., Armstrong, D., Balden, M., Baluc, N., Barthe, M.F., et al., 2013\.  Recent progress in research on tungsten materials for nuclear fusion applications in europe.  Journal of Nuclear Materials 432, 482–500.
  * Rodney et al. (2003) Rodney, D., Le Bouar, Y., Finel, A., 2003.  Phase field methods and dislocations.  Acta Materialia 51, 17–30.
  * Romero et al. (2008) Romero, I., Segurado, J., LLorca, J., 2008.  Dislocation dynamics in non-convex domains using finite elements with embedded discontinuities.  Modelling and Simulation in Materials Science and Engineering 16, 035008.
  * Ruestes et al. (2014) Ruestes, C.J., Bringa, E.M., Stukowski, A., Nieva, J.R., Tang, Y., Meyers, M., 2014\.  Plastic deformation of a porous bcc metal containing nanometer sized voids.  Computational Materials Science 88, 92–102.
  * Ryu et al. (2015) Ryu, I., Cai, W., Nix, W.D., Gao, H., 2015.  Stochastic behaviors in plastic deformation of face-centered cubic micropillars governed by surface nucleation and truncated source operation.  Acta Materialia 95, 176–183.
  * Ryu et al. (2020) Ryu, I., Gravell, J., Cai, W., Nix, W.D., Gao, H., 2020\.  Intrinsic size dependent plasticity in bcc micro-pillars under uniaxial tension and pure torsion.  Extreme Mechanics Letters 40, 100901\.
  * Ryu et al. (2013) Ryu, I., Nix, W.D., Cai, W., 2013.  Plasticity of bcc micropillars controlled by competition between dislocation multiplication and depletion.  Acta Materialia 61, 3233–3241.
  * Ryu et al. (2011) Ryu, S., Kang, K., Cai, W., 2011.  Entropic effect on the rate of dislocation nucleation.  Proceedings of the National Academy of Sciences 108, 5174–5178.
  * Schneider et al. (2009) Schneider, A., Kaufmann, D., Clark, B., Frick, C., Gruber, P., Mönig, R., Kraft, O., Arzt, E., 2009\.  Correlation between critical temperature and strength of small-scale bcc pillars.  Physical Review Letters 103, 105501\.
  * Shan (2012) Shan, Z., 2012.  In situ tem investigation of the mechanical behavior of micronanoscaled metal pillars.  JOM 64, 1229–1234.
  * Shan et al. (2008) Shan, Z., Mishra, R.K., Syed Asif, S., Warren, O.L., Minor, A.M., 2008.  Mechanical annealing and source-limited deformation in submicrometre-diameter ni crystals.  Nature materials 7, 115–119.
  * Shao et al. (2014) Shao, S., Abdolrahim, N., Bahr, D.F., Lin, G., Zbib, H.M., 2014.  Stochastic effects in plasticity in small volumes.  International Journal of Plasticity 52, 117–132.
  * Shen and Spearot (2021) Shen, Y., Spearot, D.E., 2021\.  Mobility of dislocations in fenicrcocu high entropy alloys.  Modelling and Simulation in Materials Science and Engineering 29, 085017.
  * Srivastava et al. (2013) Srivastava, K., Gröger, R., Weygand, D., Gumbsch, P., 2013\.  Dislocation motion in tungsten: atomistic input to discrete dislocation simulations.  International Journal of Plasticity 47, 126–142.
  * Srivastava et al. (2021) Srivastava, P., Jiang, K., Cui, Y., Olivera, E., Ghoniem, N., Gupta, V., 2021\.  The influence of nano/micro sample size on the strain-rate sensitivity of plastic flow in tungsten.  International Journal of Plasticity 136, 102854.
  * Stukowski et al. (2015) Stukowski, A., Cereceda, D., Swinburne, T.D., Marian, J., 2015\.  Thermally-activated non-schmid glide of screw dislocations in w using atomistically-informed kinetic monte carlo simulations.  International Journal of Plasticity 65, 108–130.
  * Tang et al. (2008) Tang, H., Schwarz, K., Espinosa, H., 2008.  Dislocation-source shutdown and the plastic behavior of single-crystal micropillars.  Physical Review Letters 100, 185503\.
  * Uchic and Dimiduk (2005) Uchic, M.D., Dimiduk, D.M., 2005\.  A methodology to investigate size scale effects in crystalline plasticity using uniaxial compression testing.  Materials Science and Engineering: A 400, 268–278.
  * Uchic et al. (2004) Uchic, M.D., Dimiduk, D.M., Florando, J.N., Nix, W.D., 2004\.  Sample dimensions influence strength and crystal plasticity.  Science 305, 986–989.
  * Vondřejc et al. (2014) Vondřejc, J., Zeman, J., Marek, I., 2014.  An fft-based galerkin method for homogenization of periodic media.  Computers & Mathematics with Applications 68, 156–173.
  * Wang et al. (2022) Wang, C., Sun, C., Cai, W., Qian, L., Guo, X., Fu, M., 2022.  Evolution of partial dislocation slip–mediated deformation twins in single crystals: a discrete dislocation plasticity model and an analytical approach.  International Journal of Plasticity 152, 103230.
  * Wang et al. (2001) Wang, Y.U., Jin, Y., Cuitino, A., Khachaturyan, A., 2001\.  Nanoscale phase field microelasticity theory of dislocations: model and 3d simulations.  Acta Materialia 49, 1847–1857.
  * Weinberger and Cai (2008) Weinberger, C.R., Cai, W., 2008\.  Surface-controlled dislocation multiplication in metal micropillars.  Proceedings of the National Academy of Sciences 105, 14304–14307.
  * Wijnen et al. (2021) Wijnen, J., Peerlings, R.H., Hoefnagels, J.P., Geers, M.G., 2021\.  A discrete slip plane model for simulating heterogeneous plastic deformation in single crystals.  International Journal of Solids and Structures 228, 111094.
  * Willot (2015) Willot, F., 2015.  Fourier-based schemes for computing the mechanical response of composites with accurate local fields.  Comptes Rendus Mécanique 343, 232–245.
  * Yaghoobi and Voyiadjis (2016) Yaghoobi, M., Voyiadjis, G.Z., 2016\.  Size effects in fcc crystals during the high rate compression test.  Acta Materialia 121, 190–201.
  * Yang et al. (2021) Yang, W., Ruestes, C.J., Li, Z., Abad, O.T., Langdon, T.G., Heiland, B., Koch, M., Arzt, E., Meyers, M.A., 2021.  Micro-mechanical response of ultrafine grain and nanocrystalline tantalum.  journal of materials research and technology 12, 1804–1815.
  * Yu et al. (2021) Yu, Q., Martínez, E., Segurado, J., Marian, J., 2021\.  A stochastic solver based on the residence time algorithm for crystal plasticity models.  Computational Mechanics 68, 1369–1384.
  * Zeman et al. (2017) Zeman, J., de Geus, T.W., Vondřejc, J., Peerlings, R.H., Geers, M.G., 2017.  A finite element perspective on nonlinear fft-based micromechanical simulations.  International Journal for Numerical Methods in Engineering 111, 903–926.
  * Zhang et al. (2006) Zhang, H., Schuster, B.E., Wei, Q., Ramesh, K.T., 2006\.  The design of accurate micro-compression experiments.  Scripta Materialia 54, 181–186.
  * Zhang et al. (2012) Zhang, L., Ohmura, T., Sekido, K., Hara, T., Nakajima, K., Tsuzaki, K., 2012\.  Dislocation character transition and related mechanical response in a body-centered cubic single crystal.  Scripta Materialia 67, 388–391.
  * Zhang and Aifantis (2011) Zhang, X., Aifantis, K., 2011\.  Interpreting strain bursts and size effects in micropillars using gradient plasticity.  Materials Science and Engineering: A 528, 5036–5043.
  * Zhu et al. (2008) Zhu, T., Li, J., Samanta, A., Leach, A., Gall, K., 2008\.  Temperature and strain-rate dependence of surface dislocation nucleation.  Physical Review Letters 100, 025502\.