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

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.


I. INTRODUCTION
Due to its thermal and electronic transport properties, it is a promising approach to use nanostructured silicon as thermoelectric material. 1 The efficiency of a thermoelectric device at temperature T increases with the material's figure of merit zT , where α is the Seebeck coefficient, σ the electrical conductivity, and κ 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 nm and longer. 2 Electrons, on the other hand, have effective mean free paths on a much shorter scale. 1 Due to advances in processing, highly doped polycrystalline silicon with crystallite sizes in the order of 25 nm can be prepared from nanoparticles. 3 For structures on the nano scale, the phonon scattering is expected to change, making phonon-surface scattering the most important mechanism. 1 The expected significant reduction in thermal conductivity was achieved by producing silicon materials consisting of nanosized crystalline grains. 3 While the thermal conductance, and thus the average influence of grain boundaries (GBs) in polycrystalline silicon films, was measured experimentally, 4 the investigation of actual individual grain boundaries is challenging. Therefore, molecular dynamics (MD) simulations can be a) Electronic mail: kevin.schroeer@uni-due.de of great help to elucidate transport processes and estimate the magnitude and behavior of thermal material properties. [5][6][7][8][9][10] 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, 7,9,11 we expect interfacial resistances in the order of a few K m 2 /GW and an increase with twist angle up to at least 40 • . 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, 12 which is well characterized and was shown to produce realistic results in previous cases. 7,10,11,13 For our simulations, we used the well established program LAMMPS. 14 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 28 Si. Particlewall 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 H ) at one end is heated by adding a certain amount of kinetic energy ∆K per time step ∆t to the atoms by rescaling their velocity vectors accordingly. The slab at the other end (thickness L C = L 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 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 int of the contained inner degrees of freedom N int at a given time step, where k 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 K equilibrium positions of the silicon lattice and are given random velocities corresponding to a Maxwell-Boltzmann distribution with a temperature of 300 K. Thereafter, the MD simulation with a time step of 1 fs consists of two phases: First, the system is equilibrated for 3 ns at 300 K using the Nosé-Hoover-thermostat 15 with constant volume V and constant number of particles N (N V T -ensemble). For large systems (L > 170 a 0 or A ≥ 200 a 2 0 ) this equilibration time is extended to 4 ns. Second, the heat flux is established and maintained for 6 ns. In this second phase, the system is given a time of 3 ns to reach the stationary state and the last 3 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 GW/m 2 ).

A. 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 κ is determined from Fourier's law (in one dimension): If Fourier's law is valid and j is constant, we expect a linear temperature profile, provided κ 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, 6,10 and accounted for by the modified phonon scattering at the heating and cooling regions.
Therefore, the profiles are fitted in the linear range and κ 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 rep-  resentation, 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.

B. 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 ϑ, 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). 16 The former are called Σ grain boundaries. Σ 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. 7 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, 17 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 K . In the presence of a heat flux, it leads to a temperature discontinuity ∆T at the interface, 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 A), 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.

A. Finite size effects on the thermal conductivity
The mono-crystalline, cylindrical system is simulated for four different wire lengths from 23.9 nm to 191 nm with a constant diameter of 6.13 nm. Figure 7 shows that the thermal conductivity increases with length and (for our system) lies between 11.6 W m −1 K −1 and 40.9 W m −1 K −1 . Schelling et al. reported a similar behavior (increase of κ 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 nm and 12.26 nm for a constant length of 47.8 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 W m −1 K −1 and 22.2 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 nanocrystalline media. 6,10,18 By comparing our results to experimental values of κ ≈ 250 W m −1 K −1 for bulk 28 Si 19 , we can confirm a reduction of the thermal conductivity in nanowires by one to two orders of magnitude.
Da Cruz et al. 10 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. Comparing the cuboidal setup with our cylindrical system of same length and cross sectional area, i.e. L = 23.9 nm and D = 6.13 nm, we note a reduction in the thermal conductivity to 11.6 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 ∼ 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. 10 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. 18 investigated isotopically natural, intrinsic, crystalline Si nanowires with diameters of 22 nm and a length of several microns, measuring thermal conductivities of about 7 W m −1 K −1 at 300 K. Moreover, the results of Boukai et al. 20 indicated experimental values of about 0.8 W m −1 K −1 for highly doped Si nanowires with cross sectional areas of 10 nm × 20 nm, and of roughly 3.5 W m −1 K −1 for A = 20 nm × 20 nm. Here, the wires were also several microns long.
We compare those values with κ = 22.2 W m −1 K −1 for our largest simulated diameter of 12.3 nm (and L = 47.8 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 K. 13 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 C).

B. 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 C, twist angle variations on the thermal boundary resistance. The 36.87 • Σ5 GB system with D = 6.13 nm is simulated for five different lengths from 23.9 nm to 382 nm, leading to a Kapitza resistance in the range from 1.46 K m 2 /GW to 1.53 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. 7 used NEMD with periodic boundary conditions to study heat conduction in Si Σ grain boundaries. They also report very little effect of length variations on the Kapitza resistance for the 43.6 • (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 • and a constant length of 47.8 nm, we use five different diameters from 4.33 nm to 17.3 nm. As can be seen in Fig. 10, the Kapitza resistance shows a decreasing and converging trend with increasing cross sectional area. Plotting R 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 with A = π D 2 /4, and a coefficient γ = 7.26 nm 2 m 2 K/GW, by the high quality of the fit.
Furthermore, we extrapolate to infinite cross sectional area (R inf = 1.23 ± 0.02 K m 2 /GW) to investigate the influence of the fixed boundary conditions. Our argument is that for A → ∞, 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 I shows that the deviation of the Kapitza resistance from R inf is 5% and less for diameters of D > 12 nm. We propose that for larger diameters, the differences in R 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 nm and L = 47.8 nm to investigate the twist angle's influence in Sec. III C. 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 • GB. Nevertheless, we expect similar behavior for other twist angles, and use the same scaling law as a first approach.

C. Dependence of Kapitza resistance on the grain boundary angle
We investigate 25 twist angles from 3.75 • to 93.75 • in equal steps of 3.75 • . Due to the twofold symmetry of the diamond lattice, angles of 90 • ± φ have equivalent configurations and should lead to the same results. This is verified for our simulations by comparing the 86.25 • and 93.75 • angle systems, which indeed yield the same Kapitza resistance of 1.48 W m −1 K −1 (see Fig. 11). Additionally, we simulate eight GB systems with coincidence site lattices. Here, the five smallest Σ values are chosen for investigation. Thus, including the 90 • Σ1 GB, we study nine Σ grain boundaries (see Table II) and a total of 33 GB systems. The determined values of the thermal boundary resistance for all angles are summarized in Table III.  Within the error bars, Fig. 11 shows a monotonous increase of the Kapitza resistance with twist angle up to 41.25 • . For larger angles, i.e., from 41.25 • to 82.5 • , R K is nearly constant with fluctuations about a mean of approximately 1.55 K m 2 /GW, as can also be seen in more detail in Fig. 12. We note a drop in R K of about 30% for the special Σ1 GB at 90 • . Besides this case, based on our results, we can not state with certainty that Σ grain boundaries have distinguished heat conduction properties compared to grain boundaries with similar angles. While the 51.13 • Σ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 Σ angles (36.87 • , 67.38 • , and 73.74 • ). 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. 7 used a method similar to ours, but with periodic boundary conditions, to simulate (001) Si twist grain boundaries for angles of 11.42 • , and 43.60 • . Their results indicate Kapitza resistances of 0.61 K m 2 /GW, and 1.25 K m 2 /GW, respectively, for T = 500 K. From theory and further simulations, the thermal boundary resistance is expected to increase for lower temperatures. 8,21 In an estimation, based on results regarding temperature dependency by Aubry et al. 21 , we assume that the values of R K should increase by about 5% to 10%, when the temperature decreases from 500 K to 300 K. Our closest angles for comparison are 11.25 • , and 45 • . Based on the finite size analysis in Sec. III B, we furthermore reduce our R K values by 17%, which leads to 0.76 K m 2 /GW, and 1.29 K m 2 /GW, respectively. Under the preceding considerations, our values agree well with those of Schelling et al. 7 , stated above.
Ju and Liang 11 studied (001) Si twist grain boundaries by NEMD simulation with periodic boundary conditions for different temperatures. Their results for T = 300 K indicate Kapitza resistances of about 1.5 K m 2 /GW, and 2.2 K m 2 /GW for angles of 16.26 • , and 36.87 • , respectively. If we extrapolate our results for these angles from Table III to infinite cross section, our values for R 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 N V T -thermostat for our fixed boundary conditions, they use a N P T -thermostat 15 (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, 22 while we use velocity rescaling. Nevertheless, the values of R 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 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 ∈ [24 nm, 380 nm] and D ∈ [4 nm, 17 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 W m −1 K −1 , which agrees well with results of da Cruz et al. 10 . Increasing this system length eightfold increased κ by a factor of 3.5. The basic configuration for the grain boundary analysis had dimensions of L ≈ 48 nm, and D ≈ 6 nm, yielding κ ≈ 19 W m −1 K −1 . Compared to experimental values for crystalline bulk 28 Si, the thermal conductivity has been reduced by one order of magnitude 19 .
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 • Σ5 GB with basic dimensions, this implied an overestimation of R K by about 17%, compared to infinite cross section. The Kapitza resistance converged quite quickly with diameter, i.e., already for D ≈ 17 nm (approx. 31 a 0 ) the deviation from the extrapolation to infinite cross section was reduced to 3%. Increasing D further to 28 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 • to 94 • (see Fig. 11). Up to 41 • , the Kapitza resistance increased with twist angle. For angles from 41 • to 86 • , we could not detect systematic behavior (see Fig. 12). The values of R K remained rather constant, fluctuating between 1.51 K m 2 /GW and 1.61 K m 2 /GW. Compared to this range, the thermal boundary resistance dropped sharply by 30% for the 90 • Σ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 Σ 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 bound-ary 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., 7 we are confident that our measured values are in a realistic range (cf. Sec. III C). 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.