Local lattice instability analysis
Original concept: bond-breaking and re-bonding necessarily
happens in inelastic deformation

Local Lattice Instability at Mode I Crack Tip:
Emergence of High Stiffness Zone in Si
Kisaragi Yashiro
Dept. Mechanical Engineering, Gifu University

Cleavage cracking

Slip deformation

Discuss local deformation by “local lattice instability”
in MD simulation
1Example: deformation process of Ni nanospecimen

Born’s stability criteria for perfect lattice
Discuss crystal stability by the positivity ofHessian, explicitly counting the interatomic potential
assuming ideal uniform lattices in continuum bodyMovie: deformation process of four surface from perfect
crystal to disordered one by slip

∂ 2U
∂ai ∂a j∂ 2U∂ai ∂a j

left: atom motion
right: colored by “local lattice instability”,
blue is stable, red unstable> 0 Stable{ ai }< 0 Unstable{ ai }Max Born (1882-1970)

2

Instability analysis by elastic stiffness coefficients

3

Stress, elastic coefficients, elastic stiffness coefficients
Elastic stiffness can be evaluated with the stress and elastic coefficientsStress and elastic coefficients are defined as
1st and 2nd order derivatives of internal energy

Stability criteria for whole crystal(Wang 1993)

Local stability at each atom point(Yashiro 2001)

45

Atomic stress and elastic coefficients in Tersoff potential
In the Tersoff potential, energy contribution of atom α is expressed as

Atomic elastic stiffness and eigenvalue analysis
Taking the Voigt symmetry into account, we represents the atomic elastic
stiffness with 6x6 matrix
instead of 4th tensor

Thus we can define atomic stress and atomic elastic coefficients as

Negative eigenvalue implies the existence of unstable deformation pathWe evaluate all the eigenvalue and eigenvector at each atom pointby using LAPACK math library

On the derivatives for three body bond order term, please see
Yashiro, Com. Mat. Sci. 112 (2016) 120–127(for 1988 Tersoff potential, however, 1989 Tersoff is rather simple so that
I would like to recommend the later version)
6Tensile simulation of center crack in Si

Simulation conditions
Potential function: Tersoff
Periodic boundary condition
Tempterature: 300K
Tensile speed:

7

Stress strain curves and movie of crack propagation

crack width: 0.3Lx
10000fs initial relaxation
Number of atoms: about 120,000
Red atom: negative 1st eigenvalue
89

Magnified view of the crack tip

Results in bcc-Fe ( [001](110) crack)(a) 1x1

Change the cell size
(different crack array)(c) 2x1

Different max. strain(b) 1x2

Movies are made with position data at every 1000fs
Unstable domain is relatively small even at the crack propagation(d) 2x2

Same distribution of unstable atoms before crack propagation
1011

Minimum 1st eigenvalue, ηα (1)-1000
Stress-2000-3000

Normalized crack length,
2w/Lx-4000

0.80.60.4
0.2

Crack length
00.05

0.1
Strain, εyy

0.15

-1000
small change

small change

-3000

-3000

-4000

10.80.60.4

Crack length
00.05

0.1

0.15

Strain, εyy-[112](111) crack, 2w=0.3Lx , T=300K[001](010) crack, 2w=0.3Lx , T=300K

10.80-1000

-1000

-1000

-2000

0.2

-2000

0.1

-3000

Crack length

0.11

0.111

0.095

Strain, εyy0.05

εyy=0.12452

ηα (1)

εyy=0.09718-3000

εyy=0.10997

0.109

Crack length
0-2000

Crack length

0.60.4[001](110) crack, 2w=0.3Lx , T=300K

00-2000

-2000

0.2

[112](111) crack, 2w=0.3Lx , T=300K
0-1000-4000

1

Detail change just at the unstable propagation

Minimum eigenvalue ηα (1) in all atoms[001](110) crack, 2w=0.3Lx , T=300K
0Normalized crack length,2w/Lx

Minimum 1st eigenvalue, ηα (1)[001](010) crack, 2w=0.3Lx , T=300K
0Normalized crack length, Minimum 1st eigenvalue, ηα (1), GPa2w/Lx

Change in the minimum eigenvalue and crack length

0.10.105

Crack length0.124

0.125

0.126

Strain, εyyStrain, εyy0.15

The eigenvalue and crack length are evaluated at every 10fs position
data (Δεyy=0.00001)

Strain, εyy

The minimum eigenvalue is selected from all the atom for each time step
The eigenvalue and crack length are evaluated at every 1000fs position
data (Δεyy=0.001)¾ large negative eigenvalue emerges before crack propagation¾ surface “negative” eigenvalues are almost zero

12

Position of large negative eigenvalue and eigenvector

13High stable area at crack tipDeterminant of 6x6 matrix of
atomic elastic stiffness

Warm color atoms have large negative eigenvalue
They are symmetrically located at crack tip(Normalized by that in perfectlattice, i.e.
in bulk)■components of eigenvector
ε1(εxx)ε2(εyy)ε3(εzz)ε4(εyz)ε5(εzx)ε6(εxy)[001](010)

-0.56

0.69

-0.12

0.41

-0.05

0.13

[001](110)

-0.59

-0.50

0.40

0.40

-0.03

0.25

[112](111)

0.30

-0.55

-0.35

-0.09

0.61

-0.33

Cold color: high stable atom
It is never seen in Fe, andmay prevent dislocation
emission and give brittle
character to Si?

The largest components suggest[001](010) crack is mode I (yy), [001](110) mode II (xx)and [112](111) mode III (zx)14

Results in the basal plane crack
Minimum 1st eigenvalue, ηα (1), GPa

Dislocations cross the cell

2a=0.1Lx

22a=0.3Lx

0.8
00.6
0.4
-50
0.2
crack length(right axis)
-100

00.01

0.02

0.03

0.04

00.05

Applied strain, εyy

1_[1010](0001), 2a=0.3Lx50

2a=0.5Lx161

ηα (1)

000.01

0.02

0.03

0.04

0.05

Applied strain, εyy

Crack opening: no unstable atom((0001) surface is stable)
Dislocation emission and crack propagation:
burst of unstable atoms1ηα (1)

0.8

0
0.6
0.4
-50
0.2
Crack length (right axis)
-100

0

0.01

0.02

00.03

Applied strain, εyy

17Normalized crack length, 2a/Lx

Stress, σyy, GPa

3

Simulation conditions
Potential function: EAM (Liu, 1996)
crack width: 0.1Lx, 0.3Lx, 0.5Lx
Periodic boundary condition
Ideal crack: virtual force shield
Tempterature: 0.1K
Number of atoms: about 120,000~130,000
Tensile speed:_[1010](0001), 2a=0.1Lx

50

Normalized crack length, 2a/Lx_[1010](0001)crack

4

Minimum 1st eigenvalue, ηα (1), GPa

Recent application: Ideal crack in Mg

15

0.4
-50
0.2
crack length(right axis)
-100

0.6
0.4
-50
0.2
crack length(right axis)
-100

00.01

0.02

0.03

0.04

1000fs00.05

Applied strain, εyy

00.02

10.8

00.6
0.4
-50
0.2

Applied strain, εyy1ηα (1)
0.8
00.6
0.4-500.2
crack length(right axis)
-100

0

0.01

0.02

0.03

0.04

00.05

Applied strain, εyy

0

-100
0.04

0.03

Dislocations cross the cell

ηα (1)

crack length(right axis)_[1010](0001), 2a=0.1Lx50

Normalized crack length, 2a/Lx

0.6

0.8
0

Dislocations cross the cell

Minimum 1st eigenvalue, ηα (1), GPa050

1ηα (1)

Minimum 1st eigenvalue, ηα (1), GPa

0.8

Minimum 1st eigenvalue, ηα (1), GPa

ηα (1)

Normalized crack length, 2a/Lx

Minimum 1st eigenvalue, ηα (1), GPa

1_[1010](0001), 2a=0.1Lx_[1010](0001), 2a=0.1Lx

50

Dislocations cross the cell

Dislocations cross the cell

Normalized crack length, 2a/Lx_[1010](0001), 2a=0.1Lx

50

Dislocation emission process (basal plane crack, 0.1Lx)

Normalized crack length, 2a/Lx

Early stage in the crack tip instability (basal plane crack, 0.1Lx)

0.05
Applied strain, εyy

Symmetrical distribution in the negative atom lines
18

Crack propagation process (basal plane crack, 0.3Lx)(a)19

Eigenvector for crack propagation(e)_[1010](0001), 2a=0.3Lx(i)
0.6(e)-50(c)

0.4(g)

0.2(h)
Crack length (right axis)

εyy

εzz

εyz

εzx

εxx

εxy

1(g)

εyy

εzz

εyz

εzx

εxy

1

0.5

0

shieldposition-0.5

0.5
shieldposition
0-0.5(h)

00.3

0.03

0.4

0.3

negative 1st eigenvalue
(i)

0.4

Normalized position, x/Lx

Normalized position, x/Lx

Applied strain, εyynegative 2nd eigenvalue

Crack instability: large normal component εxx and εzz in the 2nd eigenvector(e)20Results in the prismatic plane crack

εyy

εzz

εyz

εzx

εxy

Components of eigenvector, ε001

0.5

-0.5

00.4

Shieldposition

0.4

00.6

-0.5

0.5

0.4

-50

0.2
-100

crack length(right axis)
0.01

0.02

0.03

0.04

00.05

_
[0001](1010)crack

4

32a=0.1Lx
2

0.8
0.6
0.4

-50

0.2
Crack length(right axis)-100

0
02a=0.3Lx1

ηα (1)
00.01

0.02

0.03

0.04

0.05

Applied strain, εyy1
2a=0.5Lx

0.5

antisym.

Normalized position, x/Lx

000.010.02

0.03

0.04

0.05

Applied strain, εyy

00.4

0.6

Applied strain, εyy

10.5

0.5

0.8

010.5-0.5

0.4

η0-0.5

0.5

Dislocation motion

1α (1)

Minimum 1st eigenvalue, ηα (1), GPaεxx
10.5_[0001](1010), 2a=0.3Lx

50

Dislocation motion
Minimum 1st eigenvalue, ηα (1), GPa

0.5_[0001](1010), 2a=0.1Lx

50

1Normalized crack length, 2a/Lx

1Stress, σyy, GPa

Eigenvector for dislocation emission

2100.5-0.5

Dislocation instability: no negative 2nd eigenvalue(emerges later in the dislocation core)
0.4

-0.5

0.5

0.4

Dislocation emission rely on the unstable area?(loop like area)

0.5
2223

Normalized crack length, 2a/Lx

(f)

(a)
(b)

εxx(f)

0.8

Components of eigenvector, ε0-100(d)(d)

Normalized crack length, 2a/Lxslightlyextend!

Minimum 1st eigenvalue, ηα (1), GPa(c)1ηα (1)

Components of eigenvector, ε50(b)Results in the pyramidal plane crack

2a=0.1Lx

22a=0.3Lx
12a=0.5Lx
000.01

0.02

0.03

Applied strain, εyy

0.04

0.05

ηα (1)

0.8
0.6
-50
0.4
0.2-100
Crack length (right axis)
0Atomic elastic stiffness (AES) could be a powerful tool, suchas central symmetry parameter, to find out local deformation

100.01

Normalized crack length, 2w/Lx

Stress, σyy, GPa

3

Minimum 1st eigenvalue in all atoms, ηα (1), GPa

4

Summary__[1010](1011)crack, 2a=0.3Lx__[1010](1011)crack

The eigenvalue of the AES always shows preceding instability
before local deformation

0The normalized determinant of the AES also reveals thathighly stable zone emerges at the crack tip in Silicon.
It is never seen in metallic system.

0.02

Applied strain, εyy2425