Table Of ContentThermal resistance of grain boundaries in silicon nanowires by nonequilibrium
molecular dynamics
Jan K. Bohrer, Kevin Schr¨oer,a) Lothar Brendel, and Dietrich E. Wolf
Department of Physics and Center for Nanointegration Duisburg-Essen (CENIDE), Universita¨t Duisburg-Essen,
D-47048 Duisburg, Germany
(Dated: 26 January 2017)
The thermal boundary resistance (Kapitza resistance) of (001) twist grain boundaries in silicon nanowires
depends on the mismatch angle. This dependence is systematically investigated by means of nonequilibrium
molecular dynamics simulations. Grain boundary systems with and without coincidence site lattice are
7 compared. The Kapitza resistance increases with twist angle up to 40◦. For larger angles, it varies only little
1
around 1.56±0.05 Km2/GW, except for a drop by 30% near the 90◦ Σ1 grain boundary. Finite size effects
0
due to the fixed outer boundary conditions of the nanowire are negligible for diameters larger than 25 nm.
2
n
a I. INTRODUCTION of great help to elucidate transport processes and esti-
J
mate the magnitude and behavior of thermal material
5
Duetoitsthermalandelectronictransportproperties, properties.5–10
2
it is a promising approach to use nanostructured silicon Mainly, we will investigate the dependence of thermal
] as thermoelectric material.1 The efficiency of a thermo- resistance at twist grain boundaries on the twist angle.
ll electric device at temperature T increases with the ma- BasedonresultsofformerMDsimulationsforsiliconand
a terial’s figure of merit zT, diamond grain boundaries,7,9,11 we expect interfacial re-
h
sistancesintheorderofafewKm2/GWandanincrease
-
s α2σ withtwistangleuptoatleast40◦. Usingsystemsoffinite
e zT = T , (1)
size with fixed boundary conditions, we also aim at sim-
m κ
ulating the thermal processes between two actual grains
. where α is the Seebeck coefficient, σ the electrical con- of a given size more accurately as for periodic boundary
t
a ductivity, and κ the thermal conductivity. As z consists conditions. On the other hand, we will need to estimate
m
only of material properties, target-oriented material de- inadequacies due to the artificial fixed atomistic walls.
- signistheconsequentapproachtoimprovetheefficiency
d
of thermoelectric generators.
n
Evidently, it is desirable to find materials with a high II. METHODS
o
Seebeck coefficient, a high electrical conductivity, and a
c
[ low thermal conductivity. The approach in this work is Thermal transport in silicon is simulated by
1 to reduce the phononic contribution to heat conduction nonequilibrium molecular dynamics (NEMD) using the
v without significantly decreasing the electronic transport. Stillinger-Weber potential with original parameters,12
3 This is possible in crystalline silicon, because 50% of the which is well characterized and was shown to produce
0 thermal transport is carried by phonons with mean free realisticresultsinpreviouscases.7,10,11,13 Foroursimula-
3 paths of 500 nm and longer.2 Electrons, on the other tions,weusedthewellestablishedprogramLAMMPS.14
7 hand, have effective mean free paths on a much shorter
To obtain rotationial symmetry, we choose cylindri-
0 scale.1 Duetoadvancesinprocessing,highlydopedpoly-
cal geometry for the investigated nanowires, as shown
.
1 crystalline silicon with crystallite sizes in the order of in Fig. 1. The fixed atomistic walls have a thickness of
0 25 nm can be prepared from nanoparticles.3 For struc-
1.2 lattice constants for the cylindrical shell and one lat-
7 turesonthenanoscale,thephononscatteringisexpected
tice constant at both ends. Both the mobile part and
1
: to change, making phonon-surface scattering the most the fixed walls consist of isotopically pure 28Si. Particle-
v important mechanism.1 The expected significant reduc-
wall interactions are described by the same potential as
i tion in thermal conductivity was achieved by produc-
X particle-particle interactions.
ing silicon materials consisting of nanosized crystalline
The following method is used to establish a constant
ar grains.3 one dimensional heat flux j through the system (along
While the thermal conductance, and thus the average
z): A slab of the nanowire (thickness L ) at one end
H
influenceofgrainboundaries(GBs)inpolycrystallinesil-
is heated by adding a certain amount of kinetic energy
icon films, was measured experimentally,4 the investiga-
∆K per time step ∆t to the atoms by rescaling their
tion of actual individual grain boundaries is challenging.
velocity vectors accordingly. The slab at the other end
Therefore, molecular dynamics (MD) simulations can be
(thickness L = L ) is cooled by the same method, ex-
C H
tractingthesameamountofkineticenergypertimestep
and thus maintaining energy conservation of the whole
system. The heating and cooling regions are scaled with
a)Electronicmail: [email protected] the system length L having always a thickness of L/22.
2
L D
g g
n n
D ti mobile li mobile
a o
e o
h c
y
x
y fixed y
z
L L
H C
(a) (b) z (c) x
FIG.1. Simulatedsystemwithfixedatomisticboundaryconditions: (a)Thesimulatedmono-crystallinenanowireofcylindrical
shapecontains28Siisotopes,whichareinitiallyplacedattheirequilibriumpositionsatT =0K. (b,c)TheSi-wireiscompletely
enclosedinacylindricalshelloffixedSi-atoms. ThedimensionsofthemobilepartarediameterDandlengthL,whichincludes
theheatingandcoolinglayers(lengthsL ),wherekineticenergyisinjectedandwithdrawnrespectively. Thez-axisisparallel
H,C
to the [001]-direction of the crystal lattice.
In stationary state, the heat flux density is homoge- andaveragethetemperatureprofiles. Theheatfluxden-
neous and amounts to sity is chosen to be the same for all investigated systems
(j =21.73 GW/m2).
∆K
j = , (2)
∆tA
where A is the cross sectional area of the cylindrical A. Measuring thermal conductivity
nanowire.
The local instantaneous temperature T is defined by At first, we investigate the thermal properties of crys-
dividingthe wireintoslabs (withthickness of one lattice talline nanowires devoid of any grain boundary (GB).
constant a ) and respectively averaging the kinetic en- Knowing the heat flux and the temperature profile T(z)
0
ergy K of the contained inner degrees of freedom N fromthesimulation, thethermalconductivityκisdeter-
int int
at a given time step, mined from Fourier’s law (in one dimension):
2K dT
T := int , (3) j =−κ . (4)
N k dz
int B
wherek istheBoltzmannconstant. Thisinstantaneous If Fourier’s law is valid and j is constant, we expect
B
temperatureisthenaveragedoveratimeintervaltopro- a linear temperature profile, provided κ does not vary
duceatemperatureprofileT(z),which,inthestationary toomuchspatially. Figure2showsthatlinearityisgiven
state,doesnotdependontimeexceptforstatisticalfluc- in the middle of the nanowire. At the ends, where the
tuations. heating and cooling takes place, the system is not in lo-
Initially, the atoms are placed at the 0 K equilibrium cal equilibrium and Fourier’s law does not apply. In the
positions of the silicon lattice and are given random ve- heat source and sink, the system is externally perturbed
locities corresponding to a Maxwell-Boltzmann distribu- by rescaling the velocities and can not reach local equi-
tion with a temperature of 300 K. Thereafter, the MD librium, which also influences the adjacent regions. The
simulationwithatimestepof1fsconsistsoftwophases: nonlinear behavior was reported by other authors,6,10
First, the system is equilibrated for 3 ns at 300 K us- and accounted for by the modified phonon scattering at
ing the Nos´e-Hoover-thermostat15 with constant volume the heating and cooling regions.
V andconstantnumberofparticlesN (NVT-ensemble). Therefore,theprofilesarefittedinthelinearrangeand
For large systems (L>170 a or A≥200 a2) this equi- κ is calculated from Eq. (4). To proof consistency, we
0 0
libration time is extended to 4 ns. Second, the heat flux simulate the same system several times, only changing
is established and maintained for 6 ns. In this second the seed of the random number generator that creates
phase, the system is given a time of 3 ns to reach the the initial velocity distribution. Although the different
stationary state and the last 3 ns are used to measure profiles in Fig. 2 do not vary much in the graphical rep-
3
360
Seed #1
Seed #2
340 Seed #3
K) Seed #4 g g
e ( 320 SSeeeedd ##56 in unrotated in
atur 300 Seelidn .# f7it eat rotated ool
r
e h c
p
m 280
e
T
260 linear range
240
-25-20-15-10 -5 0 5 10 15 20 25 fixed walls
Position (nm)
FIG.3. Onehalfofthecylinder,i.e.,allatomswithpositions
FIG.2. Temperatureprofilesaveragedover3nsforaconstant
z>0,arerotatedbyanangleϑaroundan[001]-axissituated
heatfluxthroughasystemofL=47.8nmandD=6.13nm
inthecenterofthecylindercrosssection. The(001)-rotation
withoutgrainboundaryinstationarystate. Forthesamegeo-
plane lies in the x-y-plane and defines the position of the
metricalsystem,differentinitialvelocitydistributionsconsti-
grain boundary. The rotation includes the affected atoms of
tuting the same initial temperature were generated by using
the fixed walls.
differentseedsinarandomnumbergenerator. Simulationre-
sults for seven different initial velocities (seeds) show minor
variations in the temperature profiles. The thermal conduc-
tivity is calculated from the slope of the averaged profile in
the linear range. the reason for using systems with cylindrical geometry.
During the thermostat phase, an amorphous region
emerges at the grain boundary. Common neighbor anal-
ysis, as proposed by Honeycutt and Andersen,17 shows,
resentation,numericalevaluationshowsthatthethermal
which atoms are considered to be in an crystalline envi-
conductivities vary about 1%.
ronment of the diamond lattice and which have no crys-
To account for these statistical variations, we aver-
tallineorder(seeFig.4). Thethicknessoftheamorphous
age the temperature profiles of the different seeds, which
region is estimated and lies between one and three lat-
leads to a mean profile with error bars. The slope in
tice constants, depending on the twist angle and system
the linear range and its uncertainty are then determined
dimensions.
from the mean profile and its error bars.
B. Grain boundaries and Kapitza resistance
s
u
o
To implement a twist grain boundary of a certain an- h
gle, one halfofthecylindricalmono-crystallinenanowire p
r
is rotated with respect to the other, including the fixed o
m
walls, by a twist angle ϑ, as can be seen in Fig. 3. The
rotationplane(x-y)isan(001)-planeofthediamondlat- z a
tice. This prescribes the initial atom configuration. We
y
simulateGBsystemswithandwithoutacoincidencesite
lattice (CSL).16 The former are called Σ grain bound-
aries. Σ corresponds to the size of the primitive unit cell FIG. 4. After thermostating the 36.87◦ Σ5 grain boundary
of the CSL relative to the original primitive unit cell of system to 300 K, the atomic structure is examined by using
silicon. common neighbor analysis. The atoms in light gray have
diamondlatticestructure,whilethedarkgrayatomshaveno
Wechoosefixedboundaryconditions,becausewewant
crystalline order and constitute the amorphous region of the
to be able to investigate finite size effects and arbitrary
grain boundary.
grainboundarytwistangles. Forperiodicboundarycon-
ditions, only some specific angles, which create a co-
incidence site lattice, would be possible to simulate.7
Furthermore, rotation symmetry is necessary for fixed The grain boundaries constitute a resistance to ther-
boundaryconditionstoobtainequivalentregionsonboth maltransport,calledKapitzaresistanceR . Inthepres-
K
sides of the grain boundary planes. This, in the end, is enceofaheatflux,itleadstoatemperaturediscontinuity
4
∆T at the interface,
380
∆T =R j, (5) 360
K
340
which is visible in the simulated profiles (see Fig. 5). )
K
Thetemperaturejumpisdeterminedbyextrapolating e ( 320
the linear ranges on both sides of the grain boundary ur
andtakingthedifferenceofthelinearfitsattheinterface at 300
r
e
(see Fig. 5). Knowing the constant heat flux density, the p
m 280
Kapitza resistance is then calculated using Eq. (5).
e
T 260
380 240
sim.
360 lin. range lin. fit 220
lin. fit
-25-20-15-10 -5 0 5 10 15 20 25
340
K) Position (nm)
(
e 320
r
u FIG.6. Temperatureprofilesfora36.87◦ Σ5grainboundary
at 300 ∆T
r system of L = 47.8 nm and D = 6.13 nm for nine sepa-
e
p rate simulation runs with different initial velocity distribu-
m 280
tions (seeds).
e
T 260
240 lin. range
220 A. Finite size effects on the thermal conductivity
-25-20-15-10 -5 0 5 10 15 20 25
The mono-crystalline, cylindrical system is simulated
Position (nm)
for four different wire lengths from 23.9 nm to 191 nm
with a constant diameter of 6.13 nm. Figure 7 shows
FIG. 5. Measured temperature profile (sim.) averaged over
3 ns for a 36.87◦ Σ5 grain boundary system of L = 47.8 nm that the thermal conductivity increases with length
and (for our system) lies between 11.6 Wm−1K−1 and
andD=6.13nm. Thetemperaturegap∆T isdeterminedby
extrapolating the two linear ranges and taking the difference 40.9 Wm−1K−1. Schelling et al. reported a similar be-
of the linear fits at the grain boundary interface. havior (increase of κ with system length) for periodic
boundary conditions due to the reduction of the phonon
mean free path.
Again,theconsistencyoftheresultsischeckedbycom-
paring simulation runs with different initial velocity dis-
)
1
tributions. Although the graphical representation shows -K 45
no major visible differences of the temperature profiles 1 40
-m
(see Fig. 6), the numerical evaluation of the tempera-
W 35
turegapleadstosignificantdeviationsbetweentheruns.
( 30
Therefore the profiles are averaged, as for the GB-free y
systems (cf. Sec. IIA), to create a mean profile with vit 25
statistical error bars for each data point. The Kapitza ucti 20
d
resistance and its uncertainty are determined from the n 15
o
linear fits to this mean profile. c 10
al
m 5
r
e 0
h
III. RESULTS AND DISCUSSION T 0 25 50 75 100 125 150 175 200
System length (nm)
Inthissection,weanalyzefinitesizeeffectsonthether-
malpropertiesofsystemswithandwithoutgrainbound-
FIG.7. Variationofthermalconductivitywithsystemlength
ary to discuss the influence of the chosen geometry and for a system without grain boundary of constant diameter of
fixed boundary conditions. We then present our results D = 6.13 nm at T = 300 K. Error bars are smaller than
of the Kapitza resistance measurements for GB systems 1.8 % and not visible in this representation.
with 33 different twist angles.
5
Likewise, we vary the diameter six times between Forthesquarecrosssection,themobilepartofthesys-
2.45 nm and 12.26 nm for a constant length of 47.8 nm temconsistsofentirecubicunitcells. Forthecylindrical
to investigate the influence of the cross section. As geometry the unit cells near the walls will be cut into
can be seen in Fig. 8, the thermal conductivity also regions of mobile and fixed atoms. We propose that the
increases with diameter and attains values between reductioninthermalconductivityisaconsequenceofthe
15.5 Wm−1K−1 and 22.2 Wm−1K−1. However, the less regular and less periodic particle-wall bonds, which
systemlength’sinfluenceontheheatconductionismuch introduceadditionalphononscattering. Additionally, da
stronger than the diameter’s. Cruz et al.10 observed that, for simulated nanowires, the
thermal conductivity is about twice as high for fixed
) boundary conditions as for free boundary conditions in
1
-K 23 the transversal (x-y) directions.
-1m 22 Due to the challenge of characterizing nano scale sys-
W 21 tems experimentally and the imperfect material prop-
( 20 erties, comparison with real systems is difficult. Li
vity 19 et al.18 investigated isotopically natural, intrinsic, crys-
cti 18 talline Si nanowires with diameters of 22 nm and a
u length of several microns, measuring thermal conductiv-
d
n 17 ities of about 7 Wm−1K−1 at 300 K. Moreover, the re-
o
mal c 1156 saublotusto0f.B8oWukmai−e1tKa−l.120foirnhdiigchatlyeddoepxepderSimineanntaolwviraelsuewsitohf
er 14 cross sectional areas of 10 nm×20 nm, and of roughly
Th 3.5 Wm−1K−1 for A=20 nm×20 nm. Here, the wires
0 2 4 6 8 10 12 14
were also several microns long.
Diameter (nm) We compare those values with κ = 22.2 Wm−1K−1
for our largest simulated diameter of 12.3 nm (and L =
FIG. 8. Variation of thermal conductivity with diameter for 47.8 nm). Based on the preceding considerations, we
a system without grain boundary of length L = 47.8 nm at conclude an overestimation of the heat conduction in
T =300 K.
oursimulatedsystemscomparedtonaturalanddopedSi
nanowires. Reasonsaretheartificialfixedboundarycon-
ditions, the isotopic purity, the perfection of the crystal
It is well known from experiment and simulation that
lattice and material surface without impurities or oxida-
heat conduction is reduced in nano-sized and nano-
tion, and, most importantly for small systems, the dis-
crystalline media.6,10,18 By comparing our results to ex-
regard of quantum mechanical effects, which will lead to
perimentalvaluesofκ≈250Wm−1K−1 forbulk28Si19,
significant variations, even for a temperature of 300 K.13
we can confirm a reduction of the thermal conductivity
Nevertheless, NEMD can give important information
in nanowires by one to two orders of magnitude.
for thermal processes under perfect conditions and indi-
Da Cruz et al.10 used NEMD to study the thermal
cations of the thermal properties. Moreover, it is most
transportinnanowiresofsquarecrosssectionfordifferent
usefultoinvestigatetheinfluenceofparametervariations
potentials and boundary conditions. They reported sim-
bytendency,likescalingdimensions,or,inourcase,vary-
ilarfinitesizeeffects,i.e.,increaseofthethermalconduc-
ing the grain boundary twist angle (see Sec. IIIC).
tivity with length and cross sectional area A. For fixed
atomistic boundary conditions and and the Stillinger-
Weber potential, they measured κ ≈ 15.5 Wm−1K−1
for a length of 27.2 nm and cross section of 4.34 nm× B. Finite size effects on the Kapitza resistance
4.34 nm=18.8 nm2.
For a further investigation and comparison, we simu- Next, we present our results for the simulated grain
lateananowirewithsquarecrosssectionofA=29.5nm2 boundary systems, investigating the influence of length,
and L = 23.9 nm, leading to κ = 15.2 Wm−1K−1. Ex- diameter,and,inSec.IIIC,twistanglevariationsonthe
trapolating to the dimensions used by da Cruz et al. thermal boundary resistance. The 36.87◦ Σ5 GB system
(stated above), we obtain κ=15.9 Wm−1K−1, which is with D = 6.13 nm is simulated for five different lengths
in good agreement with their result of 15.5 Wm−1K−1. from 23.9 nm to 382 nm, leading to a Kapitza resistance
Comparingthecuboidalsetupwithourcylindricalsys- in the range from 1.46 Km2/GW to 1.53 Km2/GW (see
tem of same length and cross sectional area, i.e. L = Fig. 9). Although the shortest and longest system have
23.9 nm and D = 6.13 nm, we note a reduction in the highest and lowest resistance, respectively, we can
the thermal conductivity to 11.6 Wm−1K−1. Thus, by not conclude a systematic trend in the light of the given
solelychangingtheshapeofthecrosssectionfromsquare uncertainties. In fact, except for the highest value, the
tocircular,theheatconductionisreducedby∼25%. We deviationfromthemeaniscoveredbytheerrorbarsand
attribute this to inadequacies of the MD simulation, es- thehighestdeviationamountsto2.4%(seeFig.9). While
pecially in combination with fixed boundary conditions: we can not exclude a possible influence of the system
6
lengthontheKapitzaresistance, theresultsstateatany
1.8
rate that it is of minor importance for the considered lin. fit
dimensions. Schelling et al.7 used NEMD with periodic 1.7
boundary conditions to study heat conduction in Si Σ )
W
grain boundaries. They also report very little effect of G 1.6
length variations on the Kapitza resistance for the 43.6◦ 2/
m 1.5
(001) GB system. K
(K 1.4
R
GW) 1.6 mean 1.3 Rinf = (1.23 ± 0.02) K m2/GW
2m/ 1.55 1.2
0 0.01 0.02 0.03 0.04 0.05 0.06
K
e ( 1.5 D-2 (nm-2)
c
n
a
st 1.45
si FIG.10. VariationoftheKapitzaresistanceRK withdiame-
re terD fora36.87◦ Σ5grainboundarysystemwithaconstant
a 1.4 lengthofL=47.8nmatT =300K. TheKapitzaresistance
z
pit shows a linear behavior with inverse squared diameter. The
a 1.35 value for infinite cross section may be estimated by extrapo-
K
lation.
0 50 100 150 200 250 300 350 400
System length (nm)
FIG. 9. Variation of the Kapitza resistance with system TABLE I. Kapitza resistance of the simulated 36.87◦ Σ5 GB
length for a 36.87◦ Σ5 grain boundary system with a con- systems with a constant length of L = 47.8 nm for five dif-
stant diameter of D=6.13 nm at T =300 K. ferent cross sections. We include the deviation from the ex-
trapolated resistance R =1.23 Km2/GW for infinite cross
inf
section.
A Variation of the diameter, on the other hand, leads D (nm) A (nm2) R (Km2/GW) (R −R )/R
K K inf K
to a strong and systematic behavior of the Kapitza re- 4.33 14.7 1.74±0.04 0.29
sistance. For the same twist angle of 36.87◦ and a con- 6.13 29.5 1.49±0.02 0.17
stant length of 47.8 nm, we use five different diameters 8.67 59.0 1.33±0.02 0.08
from 4.33 nm to 17.3 nm. As can be seen in Fig. 10, 12.3 118 1.29±0.01 0.05
the Kapitza resistance shows a decreasing and converg- 17.3 236 1.27±0.01 0.03
ing trend with increasing cross sectional area. Plotting
R versus inverse squared diameter leads to data that is
K
ingoodagreementwithalinearfit. Althoughwedidnot However,duetothehighcomputationaleffortforlarge
establishatheorytoexplainthisbehavior,wejustifythe systems,weuseD =6.13nmandL=47.8nmtoinvesti-
approach gatethetwistangle’sinfluenceinSec.IIIC.Ifnecessary,
wearethenabletoextrapolatetoinfinitecrosssectionby
γ
R =R + , (6) reducing the boundary resistance by about 17%. The fi-
K inf A
nite size analysis was done only for the presented 36.87◦
GB. Nevertheless, we expect similar behavior for other
with A = πD2/4, and a coefficient γ =
twist angles, and use the same scaling law as a first ap-
7.26 nm2m2K/GW, by the high quality of the fit.
proach.
Furthermore, we extrapolate to infinite cross sectional
area(R =1.23±0.02Km2/GW)toinvestigatethein-
inf
fluence of the fixed boundary conditions. Our argument
is that for A → ∞, all types of lateral boundary condi- C. Dependence of Kapitza resistance on the grain
tions (free, fixed, and periodic in x-y-directions) should boundary angle
convergetothesameresults. Fromtheextrapolation,we
can estimate how big the modification due to the fixed We investigate 25 twist angles from 3.75◦ to 93.75◦ in
boundaries can be at the most. Table I shows that the equal steps of 3.75◦. Due to the twofold symmetry of
deviation of the Kapitza resistance from R is 5% and the diamond lattice, angles of 90◦ ±φ have equivalent
inf
less for diameters of D > 12 nm. We propose that for configurations and should lead to the same results. This
largerdiameters,thedifferencesinR causedbydifferent is verified for our simulations by comparing the 86.25◦
K
lateral boundary conditions should become sufficiently and 93.75◦ angle systems, which indeed yield the same
small. Kapitza resistance of 1.48 Wm−1K−1 (see Fig. 11).
7
1.8 TABLE III. Determined values of the Kapitza resistance for
the 33 investigated (001) GB systems with twist angles ϑ
at T = 300 K. Dimensions of the simulated nanowires are
1.6 L=47.8 nm and D=6.13 nm.
ϑ (◦) R (Km2/GW) ϑ (◦) R (Km2/GW)
K K
3.75 0.34±0.05 48.75 1.55±0.02
1.4 7.50 0.84±0.05 52.50 1.53±0.02
11.25 0.91±0.05 53.13 1.58±0.03
W) 15.00 1.12±0.03 56.25 1.54±0.02
G 1.2 16.26 1.16±0.02 60.00 1.55±0.05
2/ 18.75 1.17±0.03 61.93 1.59±0.06
m
22.50 1.29±0.02 63.75 1.58±0.05
K
( 1 22.62 1.28±0.02 67.38 1.52±0.06
e 26.25 1.34±0.05 67.50 1.51±0.05
c
n 28.07 1.42±0.03 71.25 1.61±0.05
a
sist 0.8 3330..7050 11..5403±±00..0022 7753..0704 11..6505±±00..0055
e
r 36.87 1.49±0.02 78.75 1.53±0.05
a
z 37.50 1.49±0.03 82.50 1.52±0.05
pit 0.6 41.25 1.55±0.02 86.25 1.48±0.06
a
K 45.00 1.55±0.02 90.00 1.08±0.04
48.75 1.55±0.02 93.75 1.48±0.05
0.4
0.2 Withintheerrorbars, Fig.11showsamonotonousin-
crease of the Kapitza resistance with twist angle up to
41.25◦. For larger angles, i.e., from 41.25◦ to 82.5◦, R
K
0
is nearly constant with fluctuations about a mean of ap-
0 10 20 30 40 50 60 70 80 90 100 proximately 1.55 Km2/GW, as can also be seen in more
GB angle (degree) detail in Fig. 12. We note a drop in RK of about 30%
forthespecialΣ1GBat90◦. Besidesthiscase,basedon
our results, we can not state with certainty that Σ grain
FIG. 11. Dependence of the Kapitza resistance on the grain
boundaries have distinguished heat conduction proper-
boundarytwistangleforasystemlengthofL=47.8nmand
ties compared to grain boundaries with similar angles.
diameter of D=6.13 nm at T =300 K.
Whilethe51.13◦ Σ5GBseemstoconstitutealocalmax-
imum in thermal boundary resistance, this could easily
be caused by statistical uncertainties or simulation in-
Additionally, we simulate eight GB systems with co- accuracies, considering the magnitude of the error bars
incidence site lattices. Here, the five smallest Σ values (see Fig. 12). The same applies to the local minima at
are chosen for investigation. Thus, including the 90◦ Σ1 some of the other Σ angles (36.87◦, 67.38◦, and 73.74◦).
GB,westudynineΣgrainboundaries(seeTableII)and Although it is reasonable to assume modified heat con-
a total of 33 GB systems. The determined values of the duction in CSLs due to the increased periodicity of the
thermal boundary resistance for all angles are summa- superlattice, we can not confirm this behavior. The ef-
rized in Table III. fect of enhanced periodicity might be suppressed by the
finite thickness of the amorphous region and the size of
the sectional area not being big enough to contain many
CSL unit cells.
TABLEII. ThenineinvestigatedΣgrainboundaries. There
Schelling et al.7 used a method similar to ours, but
aretwotwistangles(ϑ )foreveryΣvalue,symmetricabout
1,2 with periodic boundary conditions, to simulate (001)
45◦. As a trivial case, the 0◦ mono-crystalline system has
Si twist grain boundaries for angles of 11.42◦, and
no grain boundary and is listed in parenthesis for reasons of
clarity. 43.60◦. Their results indicate Kapitza resistances of
0.61 Km2/GW, and 1.25 Km2/GW, respectively, for
Σ ϑ (◦) ϑ (◦)
1 2 T = 500 K. From theory and further simulations, the
25 16.26 73.74
thermal boundary resistance is expected to increase for
13 22.62 67.38
lowertemperatures.8,21Inanestimation,basedonresults
17 28.07 61.93
regardingtemperaturedependencybyAubryet al.21,we
5 36.87 53.13
assume that the values of R should increase by about
1 (0.00) 90.00 K
5% to 10%, when the temperature decreases from 500 K
8
IV. CONCLUSION
1.7
We used NEMD to study thermal transport in crys-
) talline28Sinanowireswithandwithout(001)twistgrain
W 1.6
G boundaries. A cylindrical geometry was necessary for
2/ GB systems with fixed boundary conditions. Dimen-
m
sions were in the range of L ∈ [24 nm,380 nm] and
K
( 1.5 D ∈ [4 nm,17 nm]. The fixed boundary conditions led
e
c to significant finite size effects. The thermal conduc-
n
a tivity increased with both length and diameter, where
sist 1.4 the former had a much stronger influence (see Figs. 7
e and 8). The shortest investigated system yielded a ther-
r
a mal conductivity of about 12 Wm−1K−1, which agrees
z
pit well with results of da Cruz et al.10. Increasing this
a 1.3
K system length eightfold increased κ by a factor of 3.5.
no CSL The basic configuration for the grain boundary analysis
CSL GB had dimensions of L ≈ 48 nm, and D ≈ 6 nm, yielding
1.2 κ ≈ 19 Wm−1K−1. Compared to experimental values
20 30 40 50 60 70 80 90 for crystalline bulk 28Si, the thermal conductivity has
been reduced by one order of magnitude19.
GB angle (degree)
Concerning the thermal Kapitza boundary resistance,
we noted very little effect due to length variations (see
FIG. 12. Kapitza resistance versus grain boundary angle.
Fig. 9). In contrast, the Kapitza resistance decreased
DetailedexcerptfromFig.11foranglesfrom22.5◦ to86.25◦.
clearly and systematically with cross sectional area. Our
The eight Σ grain boundaries (CSL GB) are marked.
results indicate a linear drop with inverse cross section
(seeFig.10). Weusedthisapproachtoextrapolatetoin-
finitecrosssection,forwhichtheinfluenceoftheartificial
to 300 K. Our closest angles for comparison are 11.25◦, fixed boundary conditions should vanish. For the 36.87◦
and45◦. BasedonthefinitesizeanalysisinSec.IIIB,we Σ5 GB with basic dimensions, this implied an overesti-
furthermorereduceourRK valuesby17%,whichleadsto mation of RK by about 17%, compared to infinite cross
0.76Km2/GW,and1.29Km2/GW,respectively. Under section. The Kapitza resistance converged quite quickly
the preceding considerations, our values agree well with with diameter, i.e., already for D ≈ 17 nm (approx.
those of Schelling et al.7, stated above. 31 a0) the deviation from the extrapolation to infinite
JuandLiang11 studied(001)Sitwistgrainboundaries crosssectionwasreducedto3%. IncreasingD furtherto
by NEMD simulation with periodic boundary conditions 28 nm decreased the deviation to 1%. For future simula-
for different temperatures. Their results for T = 300 K tions with fixed boundary conditions, we suggest to use
indicateKapitzaresistancesofabout1.5Km2/GW,and diameters of at least this magnitude.
2.2 Km2/GW for angles of 16.26◦, and 36.87◦, respec- Heat conduction was simulated for 33 different (001)
tively. Ifweextrapolateourresultsfortheseanglesfrom twist grain boundaries from 4◦ to 94◦ (see Fig. 11). Up
TableIIItoinfinitecrosssection,ourvaluesforR reach to41◦,theKapitzaresistanceincreasedwithtwistangle.
inf
onlyabout60%ofJuandLiang’s. Thedeviationsmight For angles from 41◦ to 86◦, we could not detect system-
arise due to different methods used during the equilibra- atic behavior (see Fig. 12). The values of RK remained
tion of the grain boundaries prior to the heat flux phase. ratherconstant,fluctuatingbetween1.51Km2/GW and
While we use a NVT-thermostat for our fixed boundary 1.61 Km2/GW. Compared to this range, the thermal
conditions,theyuseaNPT-thermostat15 (withconstant boundary resistance dropped sharply by 30% for the 90◦
pressure P) to relax the system with periodic boundary Σ1 GB system. In the context of the calculated uncer-
conditions. Consequently, in our system, we expect both tainties,wecanneitherstatenordenydistinguishedther-
thepressuretobehigherandamodifiedstructureofthe mal properties of the other CSL grain boundaries com-
grain boundary in the vicinity of the fixed walls. How- pared to adjacent angles. It would be desirable to simu-
ever, the differences caused by different configurations lateGBsystemswithlargerdiameterformoreangles. In
shoulddiminishwithincreasingcrosssectionalarea. An- doingso,itcouldbeanalyzedwhetheraspecialbehavior
other reason for the deviation of the results might be of Σ systems is prevented by cross sections which are to
differences in the heat flux generation, as they use the small to contain enough CSL unit cells.
Mu¨ller-Plathe method,22 while we use velocity rescaling. Compared to periodic boundary conditions, we hoped
Nevertheless, the values of R are in the same order of to reproduce more accurately the thermal processes be-
K
magnitude and our simulations yields reasonable infor- tweenindividualnanoparticleswithfinitedimensions. It
mation on the grain boundary angle’s influence on ther- would be enlightening to investigate whether any major
mal boundary resistance. shortcomingsariseduetotheusedartificialfixedbound-
9
aryconditions. Thiscouldbestbedonebycomparingto 5A.Maiti,G.D.Mahan, andS.T.Pantelides,SolidStateCom-
further simulations with free boundary conditions. Be- munications102,517(1997).
cause our extrapolated values of the Kapitza resistance 6P.K.Schelling,S.R.Phillpot, andP.Keblinski,PhysicalReview
B65,144306(2002).
agree well with those of Schelling et al.,7 we are confi-
7P. K. Schelling, S. R. Phillpot, and P. Keblinski, Journal of
dent that our measured values are in a realistic range
AppliedPhysics95,6082(2004).
(cf. Sec. IIIC). At any rate, our results yield informa- 8C. Kimmer, S. Aubry, A. Skye, and P. K. Schelling, Physical
tion on the twist angle’s qualitative influence on ther- ReviewB-CondensedMatterandMaterialsPhysics75,144105
mal boundary resistance. Concerning heat conduction (2007).
9T.Watanabe,B.Ni,S.R.Phillpot,P.K.Schelling, andP.Ke-
in thermoelectric materials composed of nano particles,
blinski,JournalofAppliedPhysics102,063503(2007).
we conclude that the relative orientation of the grains 10C.A.daCruz,K.Termentzidis,P.Chantrenne, andX.Kleber,
should be of importance. Increasing the number of high JournalofAppliedPhysics110,034309(2011).
anglegrainboundariesinpolycrystallinematerialshould 11S.H.JuandX.G.Liang,JournalofAppliedPhysics113,053513
(2013).
reduce thermal conductance of the entire system.
12F. H. Stillinger and T. A. Weber, Physical Review B 31, 5262
(1985).
13S.-C. Wang, X.-G. Liang, X.-H. Xu, and T. Ohara, Journal of
ACKNOWLEDGMENTS AppliedPhysics105,014316(2009).
14S.Plimpton,JournalofComputationalPhysics117,1(1995).
15S.Nos´e,Mol.Phys.52,255(1984).
This research has been supported by the German 16A.SantoroandA.D.Mighell,ActaCrystallographicaSectionA
Research Foundation (DFG) (grant no. WO577/10-1 29,169(1973).
and WO577/6-2) within the priority program SPP 1386 17J.D.HoneycuttandH.C.Andemen,JournalofPhysicalChem-
Nanostructured Thermoelectric Materials: Theory, istry91,4950(1987).
18D. Li, Y. Wu, P. Kim, L. Shi, P. Yang, and A. Majumdar,
Model Systems and Controlled Synthesis.
AppliedPhysicsLetters83,2934(2003).
1G.Schierning,PhysicaStatusSolidi(a)15,1235(2014). 19W.S.Capinski,H.J.Maris,E.Bauser,I.Silier,M.Asen-Palmer,
2K. T. Regner, D. P. Sellan, Z. Su, C. H. Amon, A. J. H. Mc- T. Ruf, M. Cardona, and E. Gmelin, Applied Physics Letters
71,2109(1997).
Gaughey, andJ.Malen,Naturecommunications4,1640(2013).
3G. Schierning, R. Theissmann, N. Stein, N. Petermann, 20A. I. Boukai, Y. Bunimovich, J. Tahir-Kheli, J.-K. Yu, W. A.
Goddard, andJ.R.Heath,Nature451,168(2008).
A. Becker, M. Engenhorst, V. Kessler, M. Geller, A. Beckel,
H.Wiggers, andR.Schmechel,JournalofAppliedPhysics110, 21S.Aubry,C.J.Kimmer,A.Skye, andP.K.Schelling,Physical
113515(2011). ReviewB-CondensedMatterandMaterialsPhysics78,064112
4S. Uma, A. D. McConnell, M. Asheghi, K. Kurabayashi, and (2008).
K.E.Goodson,InternationalJournalofThermophysics22,605 22F. Mu¨ller-Plathe, The Journal of Chemical Physics 106, 6082
(2001). (1997).
