Thermal resistance of grain boundaries in silicon nanowires by nonequilibrium molecular dynamics

TL;DR

This study investigates how the thermal boundary resistance at grain boundaries in silicon nanowires varies with twist angle using nonequilibrium molecular dynamics simulations, revealing a dependence on mismatch angle and boundary structure.

Contribution

It provides a systematic analysis of the twist angle dependence of Kapitza resistance in silicon nanowires, including effects of boundary coincidence and finite size considerations.

Findings

Kapitza resistance increases with twist angle up to 40°

Resistance varies little (~1.56 K m^2/GW) for larger angles, except near 90°

Finite size effects are negligible for diameters >25 nm

Abstract

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 compared. The Kapitza resistance increases with twist angle up to 40{\deg}. For larger angles, it varies only little around 1.56 $\pm$ 0.05 K m^2/GW, except for a drop by 30% near the 90{\deg} {\Sigma} 1 grain boundary. Finite size effects due to the fixed outer boundary conditions of the nanowire are negligible for diameters larger than 25 nm.

Figures

Tables

Table 1. Table 1: Kapitza resistance of the simulated 36.87 ∘ 36.87\mathrm{{}^{\circ}} Σ ​ 5 Σ 5 \Sigma 5 GB systems with a constant length of L = 47.8 ​ nm 𝐿 47.8 nm L=47.8~{}\mathrm{nm} for five different cross sections. We include the deviation from the extrapolated resistance R inf = 1.23 ​ K ​ m 2 / GW subscript 𝑅 inf 1.23 K superscript m 2 GW R_{\mathrm{inf}}=1.23~{}\mathrm{K\,m^{2}/GW} for infinite cross section.

$D\left(\mathrm{nm}\right)$	$A\left({\mathrm{nm}}^2\right)$	$R_K\left(K{m}^2/\mathrm{GW}\right)$	$\left(R_K-R_{\inf }\right)/R_K$
$4.33$	$14.7$	$1.74\pm 0.04$	0.29
$6.13$	$29.5$	$1.49\pm 0.02$	0.17
$8.67$	$59.0$	$1.33\pm 0.02$	0.08
$12.3$	$118$	$1.29\pm 0.01$	0.05
$17.3$	$236$	$1.27\pm 0.01$	0.03

Table 2. Table 2: The nine investigated Σ Σ \Sigma grain boundaries. There are two twist angles ( ϑ 1 , 2 subscript italic-ϑ 1 2 \vartheta_{1,2} ) for every Σ Σ \Sigma value, symmetric about 45 ∘ 45\mathrm{{}^{\circ}} . As a trivial case, the 0 ∘ 0\mathrm{{}^{\circ}} mono-crystalline system has no grain boundary and is listed in parenthesis for reasons of clarity.

$\Sigma$	${\vartheta }_1({}_{}{}^{\circ })$	${\vartheta }_2({}_{}{}^{\circ })$
25	16.26	73.74
13	22.62	67.38
17	28.07	61.93
5	36.87	53.13
1	(0.00)	90.00

Table 3. Table 3: Determined values of the Kapitza resistance for the 33 investigated (001) GB systems with twist angles ϑ italic-ϑ \vartheta at T = 300 ​ K 𝑇 300 K T=300~{}\mathrm{K} . Dimensions of the simulated nanowires are L = 47.8 ​ nm 𝐿 47.8 nm L=47.8~{}\mathrm{nm} and D = 6.13 ​ nm 𝐷 6.13 nm D=6.13~{}\mathrm{nm} .

$\vartheta ({}_{}{}^{\circ })$	$R_K\left(K{m}^2/\mathrm{GW}\right)$	$\vartheta ({}_{}{}^{\circ })$	$R_K\left(K{m}^2/\mathrm{GW}\right)$
$3.75$	$0.34\pm 0.05$	$48.75$	$1.55\pm 0.02$
$7.50$	$0.84\pm 0.05$	$52.50$	$1.53\pm 0.02$
$11.25$	$0.91\pm 0.05$	$53.13$	$1.58\pm 0.03$
$15.00$	$1.12\pm 0.03$	$56.25$	$1.54\pm 0.02$
$16.26$	$1.16\pm 0.02$	$60.00$	$1.55\pm 0.05$
$18.75$	$1.17\pm 0.03$	$61.93$	$1.59\pm 0.06$
$22.50$	$1.29\pm 0.02$	$63.75$	$1.58\pm 0.05$
$22.62$	$1.28\pm 0.02$	$67.38$	$1.52\pm 0.06$
$26.25$	$1.34\pm 0.05$	$67.50$	$1.51\pm 0.05$
$28.07$	$1.42\pm 0.03$	$71.25$	$1.61\pm 0.05$
$30.00$	$1.43\pm 0.02$	$73.74$	$1.55\pm 0.05$
$33.75$	$1.50\pm 0.02$	$75.00$	$1.60\pm 0.05$
$36.87$	$1.49\pm 0.02$	$78.75$	$1.53\pm 0.05$
$37.50$	$1.49\pm 0.03$	$82.50$	$1.52\pm 0.05$
$41.25$	$1.55\pm 0.02$	$86.25$	$1.48\pm 0.06$
$45.00$	$1.55\pm 0.02$	$90.00$	$1.08\pm 0.04$
$48.75$	$1.55\pm 0.02$	$93.75$	$1.48\pm 0.05$

Equations

$$zT=\frac{{\alpha }^2\sigma }{\kappa }T,$$

$$j=\frac{\Delta K}{\Delta tA},$$

$$T:=\frac{2K_{\mathrm{int}}}{N_{\mathrm{int}}{k}_{\mathrm{B}}},$$

$$j=-\kappa \frac{dT}{dz}.$$

$$\Delta T=R_{\mathrm{K}}j,$$

$$R_{\mathrm{K}}=R_{\inf }+\frac{\gamma }{A},$$

Full text

I Introduction

Due to its thermal and electronic transport properties, it is a promising approach to use nanostructured silicon as thermoelectric material.Schierning (2014) The efficiency of a thermoelectric device at temperature $T$ increases with the material’s figure of merit $zT$,

[TABLE]

where $\alpha$ is the Seebeck coefficient, $\sigma$ the electrical conductivity, and $\kappa$ the thermal conductivity. As $z$ consists only of material properties, target-oriented material design is the consequent approach to improve the efficiency of thermoelectric generators.

Evidently, it is desirable to find materials with a high Seebeck coefficient, a high electrical conductivity, and a low thermal conductivity. The approach in this work is to reduce the phononic contribution to heat conduction without significantly decreasing the electronic transport. This is possible in crystalline silicon, because $50\%$ of the thermal transport is carried by phonons with mean free paths of $500~{}\mathrm{nm}$ and longer.Regner et al. (2013) Electrons, on the other hand, have effective mean free paths on a much shorter scale.Schierning (2014) Due to advances in processing, highly doped polycrystalline silicon with crystallite sizes in the order of $25~{}\mathrm{nm}$ can be prepared from nanoparticles.Schierning et al. (2011) For structures on the nano scale, the phonon scattering is expected to change, making phonon-surface scattering the most important mechanism.Schierning (2014) The expected significant reduction in thermal conductivity was achieved by producing silicon materials consisting of nanosized crystalline grains.Schierning et al. (2011)

While the thermal conductance, and thus the average influence of grain boundaries (GBs) in polycrystalline silicon films, was measured experimentally,Uma et al. (2001) the investigation of actual individual grain boundaries is challenging. Therefore, molecular dynamics (MD) simulations can be of great help to elucidate transport processes and estimate the magnitude and behavior of thermal material properties.Maiti, Mahan, and Pantelides (1997); Schelling, Phillpot, and Keblinski (2002, 2004); Kimmer et al. (2007); Watanabe et al. (2007); da Cruz et al. (2011)

Mainly, we will investigate the dependence of thermal resistance at twist grain boundaries on the twist angle. Based on results of former MD simulations for silicon and diamond grain boundaries,Schelling, Phillpot, and Keblinski (2004); Watanabe et al. (2007); Ju and Liang (2013) we expect interfacial resistances in the order of a few $\mathrm{K\,m^{2}/GW}$ and an increase with twist angle up to at least $40\mathrm{{}^{\circ}}$. Using systems of finite size with fixed boundary conditions, we also aim at simulating the thermal processes between two actual grains of a given size more accurately as for periodic boundary conditions. On the other hand, we will need to estimate inadequacies due to the artificial fixed atomistic walls.

II Methods

Thermal transport in silicon is simulated by nonequilibrium molecular dynamics (NEMD) using the Stillinger-Weber potential with original parameters,Stillinger and Weber (1985) which is well characterized and was shown to produce realistic results in previous cases. Schelling, Phillpot, and Keblinski (2004); Wang et al. (2009); da Cruz et al. (2011); Ju and Liang (2013) For our simulations, we used the well established program LAMMPS.Plimpton (1995)

To obtain rotationial symmetry, we choose cylindrical geometry for the investigated nanowires, as shown in Fig. 1. The fixed atomistic walls have a thickness of $1.2$ lattice constants for the cylindrical shell and one lattice constant at both ends. Both the mobile part and the fixed walls consist of isotopically pure $\mathrm{{}^{28}Si}$. Particle-wall interactions are described by the same potential as particle-particle interactions.

The following method is used to establish a constant one dimensional heat flux $j$ through the system (along $z$): A slab of the nanowire (thickness $L_{\mathrm{H}}$) at one end is heated by adding a certain amount of kinetic energy $\Delta K$ per time step $\Delta t$ to the atoms by rescaling their velocity vectors accordingly. The slab at the other end (thickness $L_{\mathrm{C}}=L_{\mathrm{H}}$) is cooled by the same method, extracting the same amount of kinetic energy per time step and thus maintaining energy conservation of the whole system. The heating and cooling regions are scaled with the system length $L$ having always a thickness of $L/22$.

In stationary state, the heat flux density is homogeneous and amounts to

[TABLE]

where $A$ is the cross sectional area of the cylindrical nanowire.

The local instantaneous temperature $T$ is defined by dividing the wire into slabs (with thickness of one lattice constant $a_{0}$) and respectively averaging the kinetic energy $K_{\mathrm{int}}$ of the contained inner degrees of freedom $N_{\mathrm{int}}$ at a given time step,

[TABLE]

where $k_{\mathrm{B}}$ is the Boltzmann constant. This instantaneous temperature is then averaged over a time interval to produce a temperature profile $T(z)$, which, in the stationary state, does not depend on time except for statistical fluctuations.

Initially, the atoms are placed at the $0~{}\mathrm{K}$ equilibrium positions of the silicon lattice and are given random velocities corresponding to a Maxwell-Boltzmann distribution with a temperature of $300~{}\mathrm{K}$. Thereafter, the MD simulation with a time step of $1~{}\mathrm{fs}$ consists of two phases: First, the system is equilibrated for $3~{}\mathrm{ns}$ at $300~{}\mathrm{K}$ using the Nosé-Hoover-thermostatNosé (1984) with constant volume $V$ and constant number of particles $N$ ($NVT$-ensemble). For large systems ($L>170~{}a_{0}$ or $A\geq 200~{}a_{0}^{2}$) this equilibration time is extended to $4~{}\mathrm{ns}$. Second, the heat flux is established and maintained for $6~{}\mathrm{ns}$. In this second phase, the system is given a time of $3~{}\mathrm{ns}$ to reach the stationary state and the last $3~{}\mathrm{ns}$ are used to measure and average the temperature profiles. The heat flux density is chosen to be the same for all investigated systems ($j=21.73~{}\mathrm{GW}/\mathrm{m}^{2}$).

II.1 Measuring thermal conductivity

At first, we investigate the thermal properties of crystalline nanowires devoid of any grain boundary (GB). Knowing the heat flux and the temperature profile $T(z)$ from the simulation, the thermal conductivity $\kappa$ is determined from Fourier’s law (in one dimension):

[TABLE]

If Fourier’s law is valid and $j$ is constant, we expect a linear temperature profile, provided $\kappa$ does not vary too much spatially. Figure 2 shows that linearity is given in the middle of the nanowire. At the ends, where the heating and cooling takes place, the system is not in local equilibrium and Fourier’s law does not apply. In the heat source and sink, the system is externally perturbed by rescaling the velocities and can not reach local equilibrium, which also influences the adjacent regions. The nonlinear behavior was reported by other authors,Schelling, Phillpot, and Keblinski (2002); da Cruz et al. (2011) and accounted for by the modified phonon scattering at the heating and cooling regions.

Therefore, the profiles are fitted in the linear range and $\kappa$ is calculated from Eq. (4). To proof consistency, we simulate the same system several times, only changing the seed of the random number generator that creates the initial velocity distribution. Although the different profiles in Fig. 2 do not vary much in the graphical representation, numerical evaluation shows that the thermal conductivities vary about 1%.

To account for these statistical variations, we average the temperature profiles of the different seeds, which leads to a mean profile with error bars. The slope in the linear range and its uncertainty are then determined from the mean profile and its error bars.

II.2 Grain boundaries and Kapitza resistance

To implement a twist grain boundary of a certain angle, one half of the cylindrical mono-crystalline nanowire is rotated with respect to the other, including the fixed walls, by a twist angle $\vartheta$, as can be seen in Fig. 3. The rotation plane (x-y) is an (001)-plane of the diamond lattice. This prescribes the initial atom configuration. We simulate GB systems with and without a coincidence site lattice (CSL).Santoro and Mighell (1973) The former are called $\Sigma$ grain boundaries. $\Sigma$ corresponds to the size of the primitive unit cell of the CSL relative to the original primitive unit cell of silicon.

We choose fixed boundary conditions, because we want to be able to investigate finite size effects and arbitrary grain boundary twist angles. For periodic boundary conditions, only some specific angles, which create a coincidence site lattice, would be possible to simulate.Schelling, Phillpot, and Keblinski (2004) Furthermore, rotation symmetry is necessary for fixed boundary conditions to obtain equivalent regions on both sides of the grain boundary planes. This, in the end, is the reason for using systems with cylindrical geometry.

During the thermostat phase, an amorphous region emerges at the grain boundary. Common neighbor analysis, as proposed by Honeycutt and Andersen,Honeycutt and Andemen (1987) shows, which atoms are considered to be in an crystalline environment of the diamond lattice and which have no crystalline order (see Fig. 4). The thickness of the amorphous region is estimated and lies between one and three lattice constants, depending on the twist angle and system dimensions.

The grain boundaries constitute a resistance to thermal transport, called Kapitza resistance $R_{\mathrm{K}}$. In the presence of a heat flux, it leads to a temperature discontinuity $\Delta T$ at the interface,

[TABLE]

which is visible in the simulated profiles (see Fig. 5).

The temperature jump is determined by extrapolating the linear ranges on both sides of the grain boundary and taking the difference of the linear fits at the interface (see Fig. 5). Knowing the constant heat flux density, the Kapitza resistance is then calculated using Eq. (5).

Again, the consistency of the results is checked by comparing simulation runs with different initial velocity distributions. Although the graphical representation shows no major visible differences of the temperature profiles (see Fig. 6), the numerical evaluation of the temperature gap leads to significant deviations between the runs. Therefore the profiles are averaged, as for the GB-free systems (cf. Sec. II.1), to create a mean profile with statistical error bars for each data point. The Kapitza resistance and its uncertainty are determined from the linear fits to this mean profile.

III Results and discussion

In this section, we analyze finite size effects on the thermal properties of systems with and without grain boundary to discuss the influence of the chosen geometry and fixed boundary conditions. We then present our results of the Kapitza resistance measurements for GB systems with 33 different twist angles.

III.1 Finite size effects on the thermal conductivity

The mono-crystalline, cylindrical system is simulated for four different wire lengths from $23.9~{}\mathrm{nm}$ to $191~{}\mathrm{nm}$ with a constant diameter of $6.13~{}\mathrm{nm}$. Figure 7 shows that the thermal conductivity increases with length and (for our system) lies between $11.6~{}\mathrm{W\,m^{-1}\,K^{-1}}$ and $40.9~{}\mathrm{W\,m^{-1}\,K^{-1}}$. Schelling et al. reported a similar behavior (increase of $\kappa$ with system length) for periodic boundary conditions due to the reduction of the phonon mean free path.

Likewise, we vary the diameter six times between $2.45~{}\mathrm{nm}$ and $12.26~{}\mathrm{nm}$ for a constant length of $47.8~{}\mathrm{nm}$ to investigate the influence of the cross section. As can be seen in Fig. 8, the thermal conductivity also increases with diameter and attains values between $15.5~{}\mathrm{W\,m^{-1}\,K^{-1}}$ and $22.2~{}\mathrm{W\,m^{-1}\,K^{-1}}$. However, the system length’s influence on the heat conduction is much stronger than the diameter’s.

It is well known from experiment and simulation that heat conduction is reduced in nano-sized and nano-crystalline media.Schelling, Phillpot, and Keblinski (2002); Li et al. (2003); da Cruz et al. (2011) By comparing our results to experimental values of $\kappa\approx 250~{}\mathrm{W\,m^{-1}\,K^{-1}}$ for bulk $\mathrm{{}^{28}Si}$ Capinski et al. (1997), we can confirm a reduction of the thermal conductivity in nanowires by one to two orders of magnitude.

Da Cruz et al.da Cruz et al. (2011) used NEMD to study the thermal transport in nanowires of square cross section for different potentials and boundary conditions. They reported similar finite size effects, i.e., increase of the thermal conductivity with length and cross sectional area $A$. For fixed atomistic boundary conditions and and the Stillinger-Weber potential, they measured $\kappa\approx 15.5~{}\mathrm{W\,m^{-1}\,K^{-1}}$ for a length of $27.2~{}\mathrm{nm}$ and cross section of $4.34~{}\mathrm{nm}\times 4.34~{}\mathrm{nm}=18.8~{}\mathrm{nm}^{2}$.

For a further investigation and comparison, we simulate a nanowire with square cross section of $A=29.5~{}\mathrm{nm}^{2}$ and $L=23.9~{}\mathrm{nm}$, leading to $\kappa=15.2~{}\mathrm{W\,m^{-1}\,K^{-1}}$. Extrapolating to the dimensions used by da Cruz et al. (stated above), we obtain $\kappa=15.9~{}\mathrm{W\,m^{-1}\,K^{-1}}$, which is in good agreement with their result of $15.5~{}\mathrm{W\,m^{-1}\,K^{-1}}$.

Comparing the cuboidal setup with our cylindrical system of same length and cross sectional area, i.e. $L=23.9~{}\mathrm{nm}$ and $D=6.13~{}\mathrm{nm}$, we note a reduction in the thermal conductivity to $11.6~{}\mathrm{W\,m^{-1}\,K^{-1}}$. Thus, by solely changing the shape of the cross section from square to circular, the heat conduction is reduced by $\sim 25\%$. We attribute this to inadequacies of the MD simulation, especially in combination with fixed boundary conditions:

For the square cross section, the mobile part of the system consists of entire cubic unit cells. For the cylindrical geometry the unit cells near the walls will be cut into regions of mobile and fixed atoms. We propose that the reduction in thermal conductivity is a consequence of the less regular and less periodic particle-wall bonds, which introduce additional phonon scattering. Additionally, da Cruz et al.da Cruz et al. (2011) observed that, for simulated nanowires, the thermal conductivity is about twice as high for fixed boundary conditions as for free boundary conditions in the transversal (x-y) directions.

Due to the challenge of characterizing nano scale systems experimentally and the imperfect material properties, comparison with real systems is difficult. Li et al.Li et al. (2003) investigated isotopically natural, intrinsic, crystalline Si nanowires with diameters of $22~{}\mathrm{nm}$ and a length of several microns, measuring thermal conductivities of about $7~{}\mathrm{W\,m^{-1}\,K^{-1}}$ at $300~{}\mathrm{K}$. Moreover, the results of Boukai et al.Boukai et al. (2008) indicated experimental values of about $0.8~{}\mathrm{W\,m^{-1}\,K^{-1}}$ for highly doped Si nanowires with cross sectional areas of $10~{}\mathrm{nm}\times 20~{}\mathrm{nm}$, and of roughly $3.5~{}\mathrm{W\,m^{-1}\,K^{-1}}$ for $A=20~{}\mathrm{nm}\times 20~{}\mathrm{nm}$. Here, the wires were also several microns long.

We compare those values with $\kappa=22.2~{}\mathrm{W\,m^{-1}\,K^{-1}}$ for our largest simulated diameter of $12.3~{}\mathrm{nm}$ (and $L=47.8~{}\mathrm{nm}$). Based on the preceding considerations, we conclude an overestimation of the heat conduction in our simulated systems compared to natural and doped Si nanowires. Reasons are the artificial fixed boundary conditions, the isotopic purity, the perfection of the crystal lattice and material surface without impurities or oxidation, and, most importantly for small systems, the disregard of quantum mechanical effects, which will lead to significant variations, even for a temperature of $300~{}\mathrm{K}$.Wang et al. (2009)

Nevertheless, NEMD can give important information for thermal processes under perfect conditions and indications of the thermal properties. Moreover, it is most useful to investigate the influence of parameter variations by tendency, like scaling dimensions, or, in our case, varying the grain boundary twist angle (see Sec. III.3).

III.2 Finite size effects on the Kapitza resistance

Next, we present our results for the simulated grain boundary systems, investigating the influence of length, diameter, and, in Sec. III.3, twist angle variations on the thermal boundary resistance. The $36.87\mathrm{{}^{\circ}}$ $\Sigma 5$ GB system with $D=6.13~{}\mathrm{nm}$ is simulated for five different lengths from $23.9~{}\mathrm{nm}$ to $382~{}\mathrm{nm}$, leading to a Kapitza resistance in the range from $1.46~{}\mathrm{K\,m^{2}/GW}$ to $1.53~{}\mathrm{K\,m^{2}/GW}$ (see Fig. 9). Although the shortest and longest system have the highest and lowest resistance, respectively, we can not conclude a systematic trend in the light of the given uncertainties. In fact, except for the highest value, the deviation from the mean is covered by the error bars and the highest deviation amounts to $2.4\%$ (see Fig. 9). While we can not exclude a possible influence of the system length on the Kapitza resistance, the results state at any rate that it is of minor importance for the considered dimensions. Schelling et al.Schelling, Phillpot, and Keblinski (2004) used NEMD with periodic boundary conditions to study heat conduction in Si $\Sigma$ grain boundaries. They also report very little effect of length variations on the Kapitza resistance for the $43.6\mathrm{{}^{\circ}}$ (001) GB system.

A Variation of the diameter, on the other hand, leads to a strong and systematic behavior of the Kapitza resistance. For the same twist angle of $36.87\mathrm{{}^{\circ}}$ and a constant length of $47.8~{}\mathrm{nm}$, we use five different diameters from $4.33~{}\mathrm{nm}$ to $17.3~{}\mathrm{nm}$. As can be seen in Fig. 10, the Kapitza resistance shows a decreasing and converging trend with increasing cross sectional area. Plotting $R_{\mathrm{K}}$ versus inverse squared diameter leads to data that is in good agreement with a linear fit. Although we did not establish a theory to explain this behavior, we justify the approach

[TABLE]

with $A=\pi\,D^{2}/4$, and a coefficient $\gamma=7.26~{}\mathrm{nm}^{2}\,\mathrm{m}^{2}\,\mathrm{K}/\mathrm{GW}$, by the high quality of the fit.

Furthermore, we extrapolate to infinite cross sectional area ($R_{\mathrm{inf}}=1.23\pm 0.02~{}\mathrm{K\,m^{2}/GW}$) to investigate the influence of the fixed boundary conditions. Our argument is that for $A\to\infty$, all types of lateral boundary conditions (free, fixed, and periodic in x-y-directions) should converge to the same results. From the extrapolation, we can estimate how big the modification due to the fixed boundaries can be at the most. Table 1 shows that the deviation of the Kapitza resistance from $R_{\mathrm{inf}}$ is $5\%$ and less for diameters of $D>12~{}\mathrm{nm}$. We propose that for larger diameters, the differences in $R_{\mathrm{K}}$ caused by different lateral boundary conditions should become sufficiently small.

However, due to the high computational effort for large systems, we use $D=6.13~{}\mathrm{nm}$ and $L=47.8~{}\mathrm{nm}$ to investigate the twist angle’s influence in Sec. III.3. If necessary, we are then able to extrapolate to infinite cross section by reducing the boundary resistance by about $17\%$. The finite size analysis was done only for the presented $36.87\mathrm{{}^{\circ}}$ GB. Nevertheless, we expect similar behavior for other twist angles, and use the same scaling law as a first approach.

III.3 Dependence of Kapitza resistance on the grain boundary angle

We investigate 25 twist angles from $3.75\mathrm{{}^{\circ}}$ to $93.75\mathrm{{}^{\circ}}$ in equal steps of $3.75\mathrm{{}^{\circ}}$. Due to the twofold symmetry of the diamond lattice, angles of $90\mathrm{{}^{\circ}}\pm\phi$ have equivalent configurations and should lead to the same results. This is verified for our simulations by comparing the $86.25\mathrm{{}^{\circ}}$ and $93.75\mathrm{{}^{\circ}}$ angle systems, which indeed yield the same Kapitza resistance of $1.48~{}\mathrm{W\,m^{-1}\,K^{-1}}$ (see Fig. 11).

Additionally, we simulate eight GB systems with coincidence site lattices. Here, the five smallest $\Sigma$ values are chosen for investigation. Thus, including the $90\mathrm{{}^{\circ}}$ $\Sigma 1$ GB, we study nine $\Sigma$ grain boundaries (see Table 2) and a total of 33 GB systems. The determined values of the thermal boundary resistance for all angles are summarized in Table 3.

Within the error bars, Fig. 11 shows a monotonous increase of the Kapitza resistance with twist angle up to $41.25\mathrm{{}^{\circ}}$. For larger angles, i.e., from $41.25\mathrm{{}^{\circ}}$ to $82.5\mathrm{{}^{\circ}}$, $R_{\mathrm{K}}$ is nearly constant with fluctuations about a mean of approximately $1.55~{}\mathrm{K\,m^{2}/GW}$, as can also be seen in more detail in Fig. 12. We note a drop in $R_{\mathrm{K}}$ of about $30\%$ for the special $\Sigma 1$ GB at $90\mathrm{{}^{\circ}}$. Besides this case, based on our results, we can not state with certainty that $\Sigma$ grain boundaries have distinguished heat conduction properties compared to grain boundaries with similar angles. While the $51.13\mathrm{{}^{\circ}}$ $\Sigma 5$ GB seems to constitute a local maximum in thermal boundary resistance, this could easily be caused by statistical uncertainties or simulation inaccuracies, considering the magnitude of the error bars (see Fig. 12). The same applies to the local minima at some of the other $\Sigma$ angles ($36.87\mathrm{{}^{\circ}}$, $67.38\mathrm{{}^{\circ}}$, and $73.74\mathrm{{}^{\circ}}$). Although it is reasonable to assume modified heat conduction in CSLs due to the increased periodicity of the superlattice, we can not confirm this behavior. The effect 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.

Schelling et al.Schelling, Phillpot, and Keblinski (2004) used a method similar to ours, but with periodic boundary conditions, to simulate (001) Si twist grain boundaries for angles of $11.42\mathrm{{}^{\circ}}$, and $43.60\mathrm{{}^{\circ}}$. Their results indicate Kapitza resistances of $0.61~{}\mathrm{K\,m^{2}/GW}$, and $1.25~{}\mathrm{K\,m^{2}/GW}$, respectively, for $T=500~{}\mathrm{K}$. From theory and further simulations, the thermal boundary resistance is expected to increase for lower temperatures.Kimmer et al. (2007); Aubry et al. (2008) In an estimation, based on results regarding temperature dependency by Aubry et al.Aubry et al. (2008), we assume that the values of $R_{\mathrm{K}}$ should increase by about $5\%$ to $10\%$, when the temperature decreases from $500~{}\mathrm{K}$ to $300~{}\mathrm{K}$. Our closest angles for comparison are $11.25\mathrm{{}^{\circ}}$, and $45\mathrm{{}^{\circ}}$. Based on the finite size analysis in Sec. III.2, we furthermore reduce our $R_{\mathrm{K}}$ values by $17\%$, which leads to $0.76~{}\mathrm{K\,m^{2}/GW}$, and $1.29~{}\mathrm{K\,m^{2}/GW}$, respectively. Under the preceding considerations, our values agree well with those of Schelling et al.Schelling, Phillpot, and Keblinski (2004), stated above.

Ju and LiangJu and Liang (2013) studied (001) Si twist grain boundaries by NEMD simulation with periodic boundary conditions for different temperatures. Their results for $T=300~{}\mathrm{K}$ indicate Kapitza resistances of about $1.5~{}\mathrm{K\,m^{2}/GW}$, and $2.2~{}\mathrm{K\,m^{2}/GW}$ for angles of $16.26\mathrm{{}^{\circ}}$, and $36.87\mathrm{{}^{\circ}}$, respectively. If we extrapolate our results for these angles from Table 3 to infinite cross section, our values for $R_{\mathrm{inf}}$ reach only about $60\%$ of Ju and Liang’s. The deviations might arise due to different methods used during the equilibration of the grain boundaries prior to the heat flux phase. While we use a $NVT$-thermostat for our fixed boundary conditions, they use a $NPT$-thermostatNosé (1984) (with constant pressure $P$) to relax the system with periodic boundary conditions. Consequently, in our system, we expect both the pressure to be higher and a modified structure of the grain boundary in the vicinity of the fixed walls. However, the differences caused by different configurations should diminish with increasing cross sectional area. Another reason for the deviation of the results might be differences in the heat flux generation, as they use the Müller-Plathe method,Müller-Plathe (1997) while we use velocity rescaling. Nevertheless, the values of $R_{\mathrm{K}}$ are in the same order of magnitude and our simulations yields reasonable information on the grain boundary angle’s influence on thermal boundary resistance.

IV Conclusion

We used NEMD to study thermal transport in crystalline $\mathrm{{}^{28}Si}$ nanowires with and without (001) twist grain boundaries. A cylindrical geometry was necessary for GB systems with fixed boundary conditions. Dimensions were in the range of $L\in[24~{}\mathrm{nm},380~{}\mathrm{nm}]$ and $D\in[4~{}\mathrm{nm},17~{}\mathrm{nm}]$. The fixed boundary conditions led to significant finite size effects. The thermal conductivity increased with both length and diameter, where the former had a much stronger influence (see Figs. 7 and 8). The shortest investigated system yielded a thermal conductivity of about $12~{}\mathrm{W\,m^{-1}\,K^{-1}}$, which agrees well with results of da Cruz et al.da Cruz et al. (2011). Increasing this system length eightfold increased $\kappa$ by a factor of $3.5$. The basic configuration for the grain boundary analysis had dimensions of $L\approx 48~{}\mathrm{nm}$, and $D\approx 6~{}\mathrm{nm}$, yielding $\kappa\approx 19~{}\mathrm{W\,m^{-1}\,K^{-1}}$. Compared to experimental values for crystalline bulk $\mathrm{{}^{28}Si}$, the thermal conductivity has been reduced by one order of magnitudeCapinski et al. (1997).

Concerning the thermal Kapitza boundary resistance, we noted very little effect due to length variations (see Fig. 9). In contrast, the Kapitza resistance decreased clearly and systematically with cross sectional area. Our results indicate a linear drop with inverse cross section (see Fig. 10). We used this approach to extrapolate to infinite cross section, for which the influence of the artificial fixed boundary conditions should vanish. For the $36.87\mathrm{{}^{\circ}}$ $\Sigma 5$ GB with basic dimensions, this implied an overestimation of $R_{\mathrm{K}}$ by about $17\%$, compared to infinite cross section. The Kapitza resistance converged quite quickly with diameter, i.e., already for $D\approx 17~{}\mathrm{nm}$ (approx. $31~{}a_{0}$) the deviation from the extrapolation to infinite cross section was reduced to $3\%$. Increasing $D$ further to $28~{}\mathrm{nm}$ decreased the deviation to $1\%$. For future simulations with fixed boundary conditions, we suggest to use diameters of at least this magnitude.

Heat conduction was simulated for 33 different (001) twist grain boundaries from $4\mathrm{{}^{\circ}}$ to $94\mathrm{{}^{\circ}}$ (see Fig. 11). Up to $41\mathrm{{}^{\circ}}$, the Kapitza resistance increased with twist angle. For angles from $41\mathrm{{}^{\circ}}$ to $86\mathrm{{}^{\circ}}$, we could not detect systematic behavior (see Fig. 12). The values of $R_{\mathrm{K}}$ remained rather constant, fluctuating between $1.51~{}\mathrm{K\,m^{2}/GW}$ and $1.61~{}\mathrm{K\,m^{2}/GW}$. Compared to this range, the thermal boundary resistance dropped sharply by $30\%$ for the $90\mathrm{{}^{\circ}}$ $\Sigma 1$ GB system. In the context of the calculated uncertainties, we can neither state nor deny distinguished thermal properties of the other CSL grain boundaries compared to adjacent angles. It would be desirable to simulate GB systems with larger diameter for more angles. In doing so, it could be analyzed whether a special behavior of $\Sigma$ systems is prevented by cross sections which are to small to contain enough CSL unit cells.

Compared to periodic boundary conditions, we hoped to reproduce more accurately the thermal processes between individual nano particles with finite dimensions. It would be enlightening to investigate whether any major shortcomings arise due to the used artificial fixed boundary conditions. This could best be done by comparing to further simulations with free boundary conditions. Because our extrapolated values of the Kapitza resistance agree well with those of Schelling et al.,Schelling, Phillpot, and Keblinski (2004) we are confident that our measured values are in a realistic range (cf. Sec. III.3). At any rate, our results yield information on the twist angle’s qualitative influence on thermal boundary resistance. Concerning heat conduction in thermoelectric materials composed of nano particles, we conclude that the relative orientation of the grains should be of importance. Increasing the number of high angle grain boundaries in polycrystalline material should reduce thermal conductance of the entire system.

Acknowledgements.

This research has been supported by the German Research Foundation (DFG) (grant no. WO577/10-1 and WO577/6-2) within the priority program SPP 1386 Nanostructured Thermoelectric Materials: Theory, Model Systems and Controlled Synthesis.