Supplementary Information
1
Simulation Methods
1.1
Molecular Dynamics simulations of bulk nanocrystalline Ni samples
We perfor
med molecular dynamics (MD) simula
tions with an embedded atom method (EAM) [
1]
potential for Ni [
2], which has been demonstrated to pr
ovide a reasonable description of the
structural an
d elastic properties of Ni [
3]. All simulations were perf
ormed by using the MD code
LAMMPS developed by Sandia National Laboratories [
4]. The time step in all simulations was
chosen to be 1 fs. Using a cube-axis oriented 3-D
periodic cell with a side
length of 1.76 nm of fcc
Ni, the elastic constants
of bulk Ni at 300 K were first calculat
ed as a reference. MD calculations
with a constant number of atoms,
constant volume, and constant pr
essure (NPT) at 0 GPa and 300 K
for 50 ps were first used to equilibrate the sample
at room temperature. To calculate the anisotropic
elastic constants,
1
1
c
and
1
2
c
, a simple uniaxial tensile strain rate of
71
2.84 10
s
was applied along
the x direction with an NVT ensemb
le at 300 K. The strains in the y-
and x- directions were both
controlled to be zero under the NVT ensemble. Using the resultant average stresses,
ij
, the elastic
constants were obtained from:
11
/
xx
c
, (S1)
12
()/
yy
zz
c
2
. (S2)
where the strain
is given by
t
, where
t
is the simulation time.
The system was then compressed in the [110] direction at th
e strain rate of
71
2.84 10
s
. As the
strain
applied in [110] is equivalent to
/2
xx
yy
,
xy
,
44
c
can be extracted as::
44
/
xy
c
. (S3)
The bulk modulus B is calculated to be:
B
1
(
c
2
c
)
[
3
11
12
To sim
ulate bulk nanocrystalline Ni
, five samples with periodic boundary conditions (PBC) in all 3
directions, with grain sizes of 2.2,
3.4, 4.6, 7 and 10.5 nm were cons
tructed. The side length of these
cubic samples ranged from 3.52 nm to 10.56 nm, such that each computational cell contained 64
grains. Energy minimization via the
conjugate gradient method was appl
ied to attain an equilibrium
configuration at 0 K [
6]. The structure was then annealed
at a tem
perature of 300 K, and the
pressure was controlled to be 0 GPa with NPT en
semble for 50 ps. A constant pressure between 0
GPa to 20 GPa was applied to the system for 200
ps at 300 K to calculat
e the isothermal bulk
5].
modulus using:
()
P
BV
V
. (S4)
Analogous steps were carried out through equilibr
ation to calculate the
Young’s modulus of bulk
nanocrystalline Ni. After relaxing the structure with an NPT ensemble at 0 GPa and 300 K for 50 ps,
a uniform uniaxial tensile loading was applied
in the x-direction by ex
tending the simulation box
length at the strain rate of
61
1.0 10
s
. During this straining, the
NPT ensemble was maintained
at 300K, and the pressure in the y and z direc
tions was kept at 0 GPa to allow the Poisson
from 0 to 1 GPa with five data points to ensure
ated from the
ading s
lope of the tensile stress-strain data.
contraction. The range of pressures was chosen to be
the sam
ple remained in the harmonic regime [
7, 8]. The Young’s modulus was calcul
lo
1.2.
Finite Elements Simulations of nanocrystalline Ni thin film samples
x
y
z
x
y
z
perform
ed. The commercial package Abaqus with
the linear hexahedral elements with reduced
Finite Elem
ents (FE) meshes for polycrystal
line film samples were generated by assembling N
x
N
x N
cubic grains, each with
2 nm dim
ensions. N
=N
=3 and N
=1, 2, 5, 10, 20 were used for the
sim
ulations. Each grain had cubi
c symmetry and was initially as
signed a random crystallographic
orientation [
9]. The elastic constants were
provided by the M
D simulations (see section 3.1). Perfect
inter-granular bonding was assumed and no e
xplicit modeling of th
e grain boundaries was
integration (C3D8R) was used for all FE simulations, As in MD simulations, the film thickness was
positioned along the z-direction. To simulate uniaxia
l stretching of an infinite film, unconstrained
on both sides, the following boundary conditions were imposed as shown in
Figure S 1:
- On face +x:
u
x
0.006 nm
u
x
0
- On face -x:
- On face +y:
- On face -y:
- On face +z:
free
- On face -z: free
u
y
uniform
u
y
0
conditions applied in x, y and z directions.
Figure S 1. An illustration of polycrystalline film
sample for the FE simulations, with boundary
These boundary conditions guaranteed that
y
z
0
(this condition is satisfied for each element
if all grains are isotropic, and as an integral on
the sample sides for the case of anisotropic grains).
he reaction force necessary to impose the
T
scribed displacement on the +x face,
F
pre
x
, was
extracted from the analysis, and th
e Young’s modulus was calculated as:
E
x
F
x
A
x
L
x
u
x
(S5)
where A
x
is the area of the +x face of the crystal and L
x
the length of the crystal along the x-
direction (
x
x
LdN
). A grain size,
2
dnm
was used in all simulations for the convenience of
imensional comparison with MD simulation results
, although the elastic pr
operties extracted by the
r su
ch a m
aterial, the isotropic
lculated value was
retical value.
1.3
Molecular Dynamics simulations of bulk and thin film
m
. The modified embedded
atomic method (MEAM) potential [
10] was used, and all
d
FE sim
ulations are obviously
size-independent.
To verify that these co
nditions yielded the correct modulus, a simple check was performed for an
isotropic material. A 3x3x1 crystal was built following the outlined procedure, but the shear
modulus c
44
was chosen as
c
44
(
c
11
c
12
)/2
to enforce isotropy. Fo
modulus predicted by Eqn (S7) (see sec. 2.1) was 136.094 GPa. The FE ca
136.093 GPa, which is essentially iden
tical to the theo
nanocrystalline Diamond samples
Unlike Ni, diamond has directional
covalent bonding with
a much higher cohesive energy of -7.37
eV/ato
sim
ulations were conducted us
ing LAMMPS code. MEAM potential
parameters are given in Table
S 1
.
Tab
Ref
le S 1.
EA
ential p
eters of di
amond. The names of parameters are described in
. [11
].
’
The M
M pot
aram
Elt ‘C
esub
7.37
lat ‘dia’
asub
1.0
z 4.0
t
(0)
1.0
ielem
ent 6
t
(1)
5.1
atwt
111
12.0
t
(2)
9.8
α
4.14
t
(3)
-0.57
b
(0)
4.1 rozero
1.0
b
(1)
3.0 ibar
0
b
(2)
1.5 rc
3.0
b
(3)
1.0 augt1 0
alat 3.567 attract(1,1) -1.0
2
Results
2.1.
Calculation of elastic constants of single
crystalline Ni and diamond
The elastic constants calculated
by MD at 300 K are given in
Table S 2. Agreement with
xperim
ent [
3, 12] was expected because the EAM functi
ons w
ere determined by fitting to the
lastic constants at 0 K.
e
e
Table S 2. Comparison of elastic constants fr
om (M)EAM calculation and Experimental data
C
11
(GPa) C
12
(GPa)
C
44
(GPa) B (GPa)
MEAM 1075.7 116.35 590.1 436.1
Diamond
Expt [
12]. 1084.4 127.0 576.6 446.1
EAM 248 148 126 181
Nickel
Expt[
3]. 247 147 125 180
The Young’s modulus along the <100> crystal direction is defined as:
E
100
11
/
11
where
ij
11
i
1
j
1
. This leads to:
E
100
1/
s
11
c
11
2
c
11
c
12
2
c
12
2
c
11
c
12
136 GPa
, where the comp
liance
tensor is
s
c
1
.
The Young’s modulus in any other
direction can be calculated as [
13]:
E
hkl
1
1
s
11
2
s
11
s
12
2
s
44
h
2
k
2
h
2
l
2
k
2
l
2
(S6)
here the unit vector defining the direction <hkl> is (h,k
,l). It is easily demonstrated that the stiffest
direction is <111>, for which the ob
tained anisotropi
c elastic modulus,
w
E
111
305 GPa
.
For a polycrystalline Ni sa
mple with no specific text
ure, it is reasonable to
assume that the modulus
an be well represented by the average of Eqn (S
6) over the full solid angl
ain
odel [
14], and noting that:
c
e. Using an isostr
(Reuss) m
h
2
k
2
h
2
l
2
k
2
l
2
1
4
(
h
2
k
2
h
2
l
2
k
2
l
2
)cos(
)
d
d
1
5
/
2
/2
(S7)
where
h
cos(
)cos(
),
k
cos(
)sin(
),
l
sin(
)
, we obtain a lower bound
for the modulus of an
otropic Ni polycrystal:
is
E
lb
3
s
11
2
s
12
s
44
5
205 GPa
(S8)
Similarly, an isostress (Voigt) model [
14] yields an upper bound:
E
ub
c
11
c
12
3
c
44
2
c
11
3
c
12
c
44
c
11
2
c
12
242 GPa
(S9)
The lower bound estimate is in good agreem
ent with typical e
xperimental data
(
E
exp
200 GPa
[15]).
2.2.
Bulk and Young’s moduli of bulk nanocrystalline Ni
The isotherm
al bulk modulus for na
nocrystalline Ni with the grain sizes of 2.2, 4.6 and 7 nm was
calculated by using Eqn (S4). The range of pressure
s was chosen to be from 0 to 1 GPa to ensure
the harmonic regime [
7, 8]. The bulk modulus was calculated
to be 165.7 GPa for the 2.2 nm-
grained sam
ple, 167.9 GPa for the 4.6 nm-grained
sample, and 179.3 GPa for the 7 nm-grained bulk
nanocrystalline Ni sample.
Figure S 2 shows the ratio of the calculated bulk modulus (normalized
by a typical experim
ental value for mi
crocrystalline bulk Ni (181 GPa [
16])) as a function of
decreasing grain size. The plot in
Figure S 2 clearly shows that bot
h the bulk and the Young’s
moduli are weakly dependent on the grain size in the range studied, with the smaller grain sizes
generally resulting in lower modu
li. For example, the 2.2 nm-grained sample has a bulk modulus
8.45% lower than the experimentally-determined one; while in th
e 7 nm-grained sample this
difference is virtually indistingu
ishable. This suggests that sma
ller-grained samples may exhibit
higher compressibility likely due
to the higher volume fractio
n of grain boundaries, where the
interatomic bonds are more compliant as compared
with the grain interior. The Young’s modulus
was also found to decrease from 174.5 GPa down to 117.1 GPa with the grain size reduction from
10.5 nm to 2.2 nm. The relative amount of Y
oung’s modulus decrease, 41.5%, was much more
substantial than the 8.45% reductio
n in the bulk modulus over the same range of grain sizes. This
might be explained by the more
compliant response of grain boundary bonds under shear loading.
The stress state in the bulk modulus simulations
was completely hydrostatic, and therefore
contained no shear components. In contrast
, the loading condition for the Young’s modulus
calculations introduced shear stress and strain com
ponents, leading to a mo
re pronounced size effect.
Grain size (nm)
0.9
MD
/E
Ex
0.8
0.7
E
0
2
4
6
8
10
0.5
0.6
1
Bulk modulus
Figure S 2. Bulk and Young’s moduli calculated
by MD simulations and normalized by their
respective experimental values vs
. grain size. Experimental data were extracted from Refs [
15, 16
].
Y
oung's modulus
Number of grains (through the thickness)
1.1
1
E
fi
/E
k
lm
bul
0
2
4
6
8
10
0.6
0.7
0.8
0.9
Grain size 2.2 nm
Figure S 3. Young’s m
odulus dependence on the film
thickness of 2.2 nm grained nanocrystalline
Ni thin film. Each data point was aver
aged from 5 samples with error bars.
2.3.
MD simulations of Young’s modulus
of n
anocrystalline Ni thin films
Figure S 4
displays the in-plane Young’s modulus nor
m
alized by the calculated Young’s modulus
for the bulk 2.2nm nanocrystalline Ni, as a functi
on of film thickness, expressed in terms of the
number of grains, for all the films simulated. Th
e Young’s moduli of nanocryst
alline Ni films with
thicknesses of 35.2 nm and below were found to be
lower than those of the
bulk nanocrystalline Ni.
Films with the thicknesses of ~10 grains had
close-to-macroscale Young’s modulus while the 1
grain-thick samples exhibited a ~
24.1% reduction in the Young’s modul
us as compared with bulk.
Number of grains (through the thickness)
E
film
bulk
/E
0.7
0.8
0.9
1
1.1
Grain size 2.2 nm
0.6
0
2
4
6
8
10
ness were perf
orm
ed as described in sec1.2 The Young’s modulus
of the thin film samples as a func
tion of the number of grains thr
ough the thickness is presented in
Figure S 5. Lower (Eqn (S8)) and upper (Eqn (S9)) analyt
ical lim
its for isotropic averages bound
the results, as expected
. Each sample has 3x3xN
z
grains, with the number
of grains through the
thickness, N
z
= 1, 2, 5, 10, 20. Each data point is the
average of 5 nominally identical simulations,
only differing by the randomly assigne
d grain orientations. Error bars
are displayed for each data
point.
Figure S 6 illustrates some representative meshes with Mises stress distributions.
Figure S 4. Young’s m
odulus dependence on the film
thickness of 2.2 nm grained nanocrystalline
Ni thin film. Each data point was aver
aged from 5 samples with error bars.
2.3.
Continuum effects (Finite El
e
ments simulation results)
To ascertain whether any continuum
effects could be partly responsible for the observed size effect
in the Young’s modulus of thin nanocrystalline Ni f
ilms, Finite Elements simulations of thin films
with 1-20 grains across the thick
Clearly evident from
Figure S 5
is the absence of any noticeable si
ze effect. This im
plies that the
strong dependence of the Young’s modulus on the
film thickness brought
to light by the MD
simulations (Figure S 4
) is not a continuum effect and is lik
ely attributed entirely
to the surface
effects described above. Somewhat
surprisingly, all FE simulations showed the elastic moduli 6-9%
larger than the isotropic average predictions. This is
likely a result of the c
onstraining conditions on
the y and z faces, which are forced to rema
in planar and parallel (although allowed to
expand/contract to guarantee an average uniaixial
stress condition), hence heavily constraining the
boundary grains in the system.
Figure S 5. Finite Elements predic
tion of the elastic modulus of the film as a function of the number
of grains through the thickness. Upper (Voigt) a
nd lower (Reuss) limits for isotropic averages
bound the results, as expected. The grain size used in
the simulation is 2nm, although the problem is
size independent. Notice the lack of discernible
size effect, confirming that
the strong dependence
of the modulus on film thic
kness from MD calculations (
Figure S 4) is not a continuum effect.
Figure S 6. Representative meshes
of polycrystalline thin films loaded in uniaxial tension. The
contours represent the stress distri
bution along the loading direction (
x
x
) in the grains. The film
thickness is in the z direction. (a) 3x3x1;
(b) 3x3x2; (c) 3x3x5;
(d) 3x3x10; (e) 3x3x20.
2.4
MD simulations of Young’s modulus of nanocrystalline Diamond
The MD simulations of tensile el
astic response of nanocrystalline
Ni revealed a strong dependence
of the Young’s modulus on both the grain size and the film thickness. The spatial profile of the local
Young’s moduli, shown in Figure 3, demonstrated
that the atomic structures within the grain
boundaries and in the vicinity of
the free surface were more comp
liant than those in the grain
interior. These findings suggest th
at the regions within a material, where the local spacing among
the atoms is larger than the equilibrium distan
ce, deform elastically more readily than the
equilibrium-positioned portions (atoms at the free surface can be regarded to be infinitely apart
from their neighbors along the outw
ard normal direction of the surf
ace). This relationship between
the atomic structure and the ensuing mechanical
response has been linke
d to the shape of the
interatomic potential curve because
the elastic properties are closely
related to the second derivative
of the potential curve with respec
t to the interatomic distance [
17]. Comparing the influence of
grain s
ize and thin film thickness on the elastic pr
operties between nanocrystalline Ni and a material,
whose bonding energetics are vastly different from
those of Ni, could help shed light on the
underlying physics behind these phenomena.
Figure S 7
shows the Young’s modulus of diamond as a f
unction of both the grai
n size and the film
thickness, calculated by both m
ethods (see secti
on 2.3). Note that comparing the stress level of
dashed lines can show th
e grain size effects on Young’s modulus
of bulk diamond. This plot reveals
that the effects of both the grain size and the film thickness also exist in nanocrystalline diamond:
Young’s modulus decreases with the reduction in either
parameter. The results
of stress-strain
method (10
-1
strain), (M1) are shown to agree w
ith those of long time NPT method (~10
-3
strain)
(M2), ensuring that the deformation in both
methods is done in an elastic regime.
Figure S 7. Young’s modulus of nanocrystalline diamond as a func
tion of grain size and film
thickness. Dashed lines repr
esent Young’s modulus of bulk forms, with the green one
corresponding to the coarse-grained polycrystalline diamond cal
culated by Reuss average in
equation (S8). In the legend, M1 st
ands for the stress-strain method
and M2 stands for the long time
NPT method.
2.3
Discussions
Young’s modulus is proportional to
the second derivative of the potent
ial energy curv
e with respect
to the interplanar spacing of crystallographic plan
es, whose normal is parall
el to the axial loading
direction. If the pote
ntial energy curve is harmonic, th
e Young’s modulus remains constant
regardless of the interplanar spacing. In real mate
rials, the shape of the potential energy curves is
usually anharmonic. The Young’s
modulus of grain boundaries,
which usually accommodate a
certain amount of free volume, and hence have hi
gher-than-equilibrium atomic
spacing, is expected
to be lower than that of the grain interior. Th
e discrepancy between the Young’s moduli in the grain
boundary vs. grain interior is corre
lated with the amount of anharm
onicity in the potential energy
curve and may explain the diffe
rent Young’s moduli reduction be
tween Ni and diamond studied
here. The finite element analysis performed in
section 3.5 does not have the information of the
potential energy curve with respect to the location in the samples. T
hus, it is not surprising that the
results of the finite element analysis
did not show any length scale effect.
To quantify the degree of harmonicity, we st
udied the potential ener
gy variation during the
separation of two rigid crystalline bl
ocks. This approach is preferred
to the more typical analysis of
two neighboring atoms, as it
better captures the Young’s modulus
of volumetric samples. We
cleaved the surfaces along the (
001), (111) and (110) cr
ystallographic planes,
and then computed
the potential energy as a functi
on of the separation distance,
(see
Figure S 8
(a)). These three
particular planes represent some of the most
common orientations, and serve as reliable model
surfaces for the potential energy quantifi
cation. The simulation block sizes were
12[100]×4[010]×10.5[001],
12[11
2]
4[
110]
10.5[111]
and
12[111]
4[11
2]
10.5[
110]
. Figure S
8(a) shows an example of such a
po
tential energy curve as a func
tion of separation distance for the
separation on the (001) plane for Ni, and
Figure S 8
(b) shows its second derivative, (
).
Both variables in the plot are normalized by their respective equilibrium values, denoted by the
subscript “0,” to avoid the effects of simulation
cell size. A greater degree of anharmonicity in the
potential curve would result in a more
pronounced drop in curvature at larger
2
2
E
pot
/
/
0
.
Figure S 8. (a) The potential energy variation of Ni
with respect to the se
paration distance between
two crystalline blocks (
) for (001) plane. (b) Curvature of the potential energy shown in (a)
as
0
22
22
/
pot
pot
EE
a function of the separation distance.
Both x- and y-axes are normalized by
their equilibrium values for
comparison between two different
materials (Ni and diamond)
Following this methodology, the second
derivatives of the potential en
ergy curves for the separation
along (001), (111), an
d (110) planes for Ni
and diamond are shown in
Figure S 9. These plots
consistently dem
onstrate that th
e curvature of the potential curve
for Ni falls off with separation
distance more rapidly than that for C. This resu
lt implies that the potential curve of Ni is more
anharmonic than that of diamond, and
that for a given
atomic spacing (
/
0
), the effects of grain
size on Young’s modulus are stronger for nanocryst
alline Ni than for nanocrystalline diamond.
These findings are consistent with the simulatio
ns performed in this work (Figure 2 (a) and
Figure S
7).
Figure S 9. Curvature of the potential energy curves as a function of normalized separation distance
along (a) (001), (b
) (111) and (c) (110) planes.
Our simulations demonstrate that
the Young’s modulus decreased wi
th decreasing film thickness,
which suggests its dependence on th
e film-thickness-to-grain size ra
tio (t/d). Similar “smaller is
weaker” trend in the yield strength has only been observed
in a size-dependent study in
nanocrystalline Ni nanopillars [
18], whereby nanocrystalline Ni–4%W
nanopillars with
grain sizes
of 60nm exhibited a “smaller is w
eaker” trend at smaller D/d’s. Recent work on the deformation of