Radial dependence of thermal transport in silicon nanowires

A radial decomposition of the heat flux in a silicon crystalline nanowire is studied with molecular dynamics (MD) and Monte Carlo (MC) simulations. Less heat flux is carried in the external layer of nanowires than in the center. The difference between the center and the surface is of the order of 50% and 30% with MD and MC simulations, respectively. As a result, a heat flux close to the surface is 30% and 15% lower than the total axial heat flux in the structure. The physical mechanism behind is analyzed from partial contribution of each atom in each phonon mode calculated from lattice dynamics. The reduction of the flux close to the surface is related to back-scattering, the amorphous-like DOS of the external layer and flattened dispersion curves, thus lower phonon group velocities. Our study points to the need for cautions analysis of experimental determination of the thermal conductivity involving contact measurements, such as scanning thermal microscopy.


Introduction
The thermal conductivity of nanowires has been a popular issue the last two decades due to the potential applications in thermal management and energy harvesting. As the fabrication methods evolve rapidly several issues emerged concerning the modification of their physical properties, in particular their thermal properties. In semiconductor nanowires, phonons are the dominant energy carriers. Several studies showed that the thermal conductivity of nanowires is not anymore an intrinsic property but it depends on the diameter, length, surface roughness, composition, doping, crystallinity and crystalline orientation, ending facets, native oxides or eventual nanoconstriction of nanowires [1]. The thermal conductivity of silicon nanowires can be reduced by a factor up to 100 depending on the above parameters. In general the severe reduction of the thermal conductivity in nanowires is related to the increase of the surface to volume fraction thus the increase of the boundary scattering, and the phonon confinement due to the reduction of the dimensionality which leads to softened phonon dispersion and reduced phonon group velocities. Donadio and Galli argued that non-propagating modes and decrease of the lifetimes of propagating modes are the reasons for the dramatic reduction of the thermal conductivity of nanowires [2].
There is a plethora of experimental measurements and simulations for the thermal conductivity of silicon nanowires, and the results are quite dispersed. Theoretical articles show a variation of the thermal conductivity from 0.1 up to 50Wm −1 K −1 depending on their diameter [2][3][4][5][6][7][8] . In most of them there is a deviation from Fourierʼs law [9]. Among the main methods for measuring the thermal conductivity of nanowires one finds the bridge technique [10][11][12] and the tip methodology [13]. While in the first method the nanowires are attached in their extremities and a cold and a hot reservoir are imposed, thus the total thermal conductivity is measured, in the second methodology a hot AFM tip approaches the nanowire in a certain position and the thermal conductivity is deduced from the contact with the external layer.
We note that in many nanostructures there is a native oxide layer at the surface which influences the heat conduction. In addition, rough surfaces should diffuse more phonons compared to smooth surfaces. The question treated here is if the surfaces conduct better or worse the heat current even if there are totally smooth Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. and free of amorphous phase shell layers including native oxides. This could have a severe influence on the thermal conductivity measurements with contact and scanning probe technique [14].
In this work we investigate the radial dependence of the heat flux in a cylindrical nanowire. All nanowires are silicon crystalline nanowires with growth direction [100] and they are free of atomic roughness or presence of amorphous shells. The article is divided in four sections. After the introduction, details of the modeling of nanowires and the simulation method are given in the second section, including the methodology to calculate the partial DOS as well as the analytical model of Dingle [15]. In the third section the results of molecular dynamics (MD) and Monte Carlo (MC) simulations are given for several nanowires with small diameters ( 11 < nm). Due to CPU constrains we studied large diameters only with the MC method. In the last section conclusions and discussions are given.

Computational details
2.1. Molecular dynamics (MD) MD simulations are performed with LAMMPS [16]. The timestep is set to 0.5fs and the Stillinger-Weber potential [17] with modified coefficients [18] is used. First, a cylindrical silicon nanowire with variable diameter d (from 5 to 11nm) and constant length L tot =240 nm is designed. Then the structure is relaxed during an equilibration stage under NVT ensemble at room temperature (T=300K) for 200ps. After this step, the system is divided in 5 concentric cylindrical layers with thickness t=d/10 parallel to the axis of the nanowire (see figure 1). The global NVT ensemble is replaced by a NVE one, and a hot and a cold thermostats are applied at both ends of the structure, with temperatures T H =330 K and T C =270 K, respectively. The axial heat flux J z in each cylindrical layer is then computed and recorded during 6ns, by means of the expression with V the volume of the considered layer, E the sum of the potential and kinetic energies, v  the atom velocity, and S the atomic stress tensor. The sum is made over all atoms contained in the layer. Every simulation is repeated several times (typically 8 times) with different initial conditions in order to reduce the uncertainty of the results.
Each thermostated region has a length of 11nm, so the distance between the hot and cold zones is L=218nm. The diameter of the nanowires studied with MD are comprised between 5 and 11nm (see table 1). With such a small characteristic dimension, the phonon mean free path (MFP) is drastically reduced and the diffusive regime is retrieved for relatively small nanowire length. Thus, we can consider that there is no significant contribution of ballistic transport between the thermostats.

MC resolution of Boltzmann Transport Equation (BTE)
The MC technique of the resolution of BTE is based on the modeling of phonon propagation and scattering from a particle point of view. Details of the algorithm can be found in our previous works [19][20][21]. The BTE is solved under the relaxation time approximation. The phonon properties (dispersions, density of states (DOS), group velocities, lifetimes) are those of bulk silicon and the lifetimes are obtained from Hollandʼs formalism [22]. The phonon properties of a nanowire usually differ from bulk ones, but it has been previously shown that it has a weak impact on thermal transport in silicon nanostructures at room temperature [23][24][25]. Moreover, it is not trivial to get the phonon properties for nanostructures, while bulk ones are widely documented in the literature. Thus, the use of bulk phonon properties in MC simulations is common and usually provides results in agreement with experimental data [5,21,26]. The surface of the cylindrical nanowire is set to purely diffuse, which is consistent with the usual roughness of experimentally elaborated nanowires (δ;1 nm) [27] and the dominant wavelength in silicon at room temperature (λ;1 nm) [28]. As in molecular dynamics simulations, a temperature gradient is applied and the cylindrical nanowire is divided in 5 concentric cylindrical layers along the nanowire axis. In MC resolution of BTE, phonons are considered as particles with a defined position. Thus, it is trivial to determine in which layer a phonon is localized at the end of a timestep. Then the heat flux can be computed every timestep in each layer summing the contribution of phonons which are in the layer at the end of the current step and normalizing by the appropriate volume.
With MC, larger diameters can be investigated (up to 20μm). With such large diameters, the distance between the thermostats has to be long enough to avoid ballistic transport by heat carriers. Thus, the nanowire length is set to different values for each studied diameter, and L varies from 1 to 10μm as d goes from 5 to 20nm (see table 1). According to the dominant phonon intrinsic MFP in bulk silicon at room temperature (∼300nm) [29], a length of 10μm is sufficient to avoid ballistic transport as the diameter is very large. The timestep is set to 1ps and the total duration of a simulation is about 1.5μs. The temperature difference is set to 2K (T H =301 K and T C =299 K).

Results and discussion
In nanoscale objects thermal conductivity is not shape independent anymore. It is obtained by dividing the total flux along the gradient direction by the cross section. Here a local thermal conductivity is also defined, from the local flux in the z-direction as a function of the radius. The radial variation of the thermal conductivity is then obtained for small diameters with both methods described in the previous section, as the ratio of the local axial heat flux to the total temperature gradient (Fourierʼs law). The results are plotted in figure 2. The same trend is observed for all diameters and with both methods: thermal conductivity is maximum at the center of the nanowire and decreases when approaching the surface, as expected from theoretical models based on electronic transport or hydrodynamic considerations [15,30,31]. From a thermal point of view, the reduction of thermal conductivity close to the surface can be explained with phonon surface scattering. At this scale, phonons propagating in the external layer are scattered more often by the nanowire surface than phonons at the center of the nanowire. Thus, thermal resistance is enhanced in the external layer, especially when backscattering is important (i.e. for diffuse phonon reflexions).
The dashed lines in figure 2 represent the usual total thermal conductivity along the nanowire axis. It is computed from the total heat flux in the structures, but it can also be retrieved with an average of the thermal conductivities of each layer, weighted by their surface. With both methods, the absolute value of the total thermal conductivity generally increases with the diameter (except for MD results with diameters of 8.1 and 10.9nm, but the difference is in the error bars). But for the same diameter, MD predicts larger thermal conductivity than MC. This could be related to the surface specularity, which is high in MD simulations [32], and zero with MC because of the diffuse boundary condition. In order to compare the results obtained with the two numerical methods for different diameters, we normalized the thermal conductivity. In this work, we focus on relative variations of thermal properties with the radius.
In figure 3, the results obtained with the two methods are gathered in order to compare them. The heat flux computed in each layer is normalized by its maximum value (the one at the center of the nanowire). As the temperature gradient is the same for all cylindrical layers, the normalized heat flux is equivalent to the normalized thermal conductivity. The trend given by both methods is the same, but MD predicts a more pronounced reduction of the thermal transport close to the surface compared to the center (∼50%) than MC (∼30%). This result is counter-intuitive, as the nanowire surface is more specular with MD, so surface scattering alone can not explain greater reduction compared to MC results. Compared to the total axial thermal transport in the nanowire, the reduction in the external layer is found to be∼30% and∼15% with MD and MC, respectively. The numerical results can be compared to an analytic solution of the BTE derived by Dingle [15] in the case of a cylindrical wire: with p the specularity parameter for the nanowire surface and Λ the dominant MFP in bulk silicon (Λ=268 nm [29]). In figure 3, the extremes cases of p=0 (fully diffuse surface) and p=1 (fully specular surface) are plotted.  When considering only specular reflections, the heat flux has no dependence on the local radius because there is no backscattering to lower thermal transport close to the surface. When the surface is considered as fully diffuse, the model exactly reproduces the MC results, as MC simulations are performed setting the specularity parameter to 0. Equation (2) gives very close curves for d=5.4 nm, d=8.1 nm and d=10.9 nm, as observed with the simulations, so only the solution for d=8.1 nm is plotted here. MD results show a greater reduction of thermal transport than the model with p=0. This is the proof that backscattering alone cannot explain the radial variation of the flux, otherwise it would suppose a negative specularity parameter.
In figure 3, we can also observe that for each method the relative reduction of the heat flux seems independent of the nanowire diameter, as predicted by the analytical model. If this reduction is really due to surface scattering, it should vanish for larger diameters, when the characteristic dimension of the structure becomes greater than the phonon intrinsic MFP, which is the same in the entire structure. In order to check this, the radial evolution of heat transport has been investigated for nanowires with large diameters with the MC method (MC 4 to MC 6 samples). The results are plotted in figure 4.
For d=1μm, the reduction of thermal transport close to the surface is still as strong as for small diameters (25%-30%). But for d=20μm, thermal transport is quasi-constant along the radial cross-section of the nanowire. A slight reduction (∼5%) is observed in the external layer, even if the diameter is very large compared to the dominant phonon intrinsic MFP. With such a diameter, the thickness of each cylindrical layer is 2μm (see table 1). In silicon at room temperature, some phonons can have MFPs of the order of the micrometer and their contribution to heat transport is non negligible [23]. Thus, in an external layer of 2 μm of a nanowire, one can still observe a reduction of the thermal conductivity due to surface scattering of these phonons. But for macroscopic dimensions, the reduction would not be visible, except very close to the surface (<10 μm).
Equation (2) is also evaluated for the large diameters with p=0 (purely diffuse nanowire surface) and compared to MC simulations in figure 4. For d=115nm and d=1μm, it is in excellent agreement with the numerical results. For d=20μm, it correctly predicts that the reduction of thermal transport vanishes when the dimension becomes very large compared to the intrinsic MFP. Doing the average of the analytical solution over the external shell for d=20μm gives a normalized heat flux of 96.1%, which is in perfect agreement with the value obtained with MC simulations (last pink diamond).
In MD, as the nanowire surface is partly specular, other effects should be taken into account to understand the strong reduction of thermal transport close to the surface. For this purpose, the partial phonon DOS and dispersion curves have been computed in each layer. The DOS is calculated by a Fourier transform of the autocorrelation of atom velocities during a simulation, a method which has been used before [33,34]. For this, we used the specden command of the Signal Processing Plugin Package [35] of Visual molecular dynamics [36]. The DOS of the different layers of a nanowire with d=5.4 nm are compared in figure 5(a).
In the internal layers, the DOS is almost exactly the same as the one of bulk crystalline silicon (the latter has been computed with the same methodology modeling a bulk system with MD). Only the layer closest to the surface has a different DOS. A great reduction of the peak associated to the transverse optical phonons is observed (at f ; 500 cm −1 ), and phonon population is globally shifted toward lower frequencies ('redshift'). Thus, in average phonons transport less energy (E w = ). Strikingly, the DOS of the external layer is very similar to the one of bulk amorphous silicon, as can be seen in figure 5(b). The thermal conductivity of amorphous silicon is very low compared to that of crystalline silicon (∼1.5Wm −1 K −1 and ∼150Wm −1 K −1 , respectively), but this is usually associated to the extremely small phonon MFP in amorphous materials (few nanometers). Nevertheless, the quasi-amorphous DOS of the external layer surely lowers thermal transport close to the surface. To analyze the radial evolution of thermal transport in silicon nanowire, we first calculated the phonon dispersions of the silicon nanowire with diameter of 5.4 nm, with the phonopy code by solving the eigenvalues of the dynamical matrix constructed from the harmonic force constants [37,38]. The phonon dispersions are depicted in figure 6. The absence of imaginary frequencies in the phonon dispersions confirms the nanowire thermodynamical stability. For all dispersion curves, there are strong coupling between acoustic and optical phonon branches at frequencies larger than 0.5 THz, leading to an intuitive expectation of strong acousticoptical phonon scattering and thus low intrinsic lattice thermal conductivity (see figure 2) compared to its bulk counterpart (κ bulk ; 150Wm −1 K −1 ) [39].
To further explore the physical mechanism behind the trend that the thermal conductivity in silicon nanowire declines strongly from center to external layer and finally reaches a minimum at its surface, the partial contribution of each atom in each phonon mode is calculated, summarized by spacial layer, and finally averaged by the number of atoms in the layer, as shown by bubbles in figure 6. This technique has been widely used to analyze the atom-phonon correlation in previous works [40][41][42]. Here, we only consider some representative contribution larger than some criteria (see caption of figure 6) but neglect the minor ones which are not critical to our analysis. Finally only the bubbles which meet this criterion are displayed (bubbles alive for each layer: N 1 = 6901, N 2 = 1639, N 3 = 152, N 4 = 50, N 5 = 2708). First, we can see from the figure that the first layer, i.e. the center one (red bubbles), contributes nearly in the whole phonon frequency range in the total phonon dispersion curves. Also, the number of bubbles is the largest one of the five, indicating that the first layer carries the main phonon vibration information of the whole nanowire and contributes most to the total thermal conductivity. From the structure aspect, the central layer is far away from the surface, therefore modes feel less the boundary conditions and thus lattice can still vibrate as in the bulk state. Second, we can observe in the same figure that the number of bubbles for the layers from the center (layer 1) to the external layer (until layer 4) decreases rapidly, e.g. the blue bubbles of layer 3 sparsely appear at the frequency range from about 2 to 10 THz, and only a few magenta ones for layer 4 are found at around f = 16 THz (flat optical phonon branches). These numbers can be seen roughly as an estimation of the intensity of lattice vibration and thus their contribution to the total thermal conductivity. This analysis agrees well with our MD and MC results. Last but not least, we notice that the number of bubbles alive for the outside surface layer (layer 5 by cyan bubbles) is quite large, but this layer thermal conductivity is the lowest of all. This can be attributed to the flat surface states in the phonon dispersion curves with almost negligible group velocity which are formed by the vibration of the surface atoms.
Finally, it is worthwhile to explain the correlation of the geometric structure of layers in the silicon nanowire and their contribution to the total phonon dispersion. The inner 4 layers contribute mainly to the bulk silicon states, although phonon dispersion curves are folded many times due to the nanowire width and the large number of atoms in a unit cell compared to the bulk state. But the surface states (from layer 5) are newly introduced to the bulk phonon dispersion curves due to the surface reconstruction, dangling bonds formed and atoms environment sudden change there. That is why we can notice almost absolute surface states at Z point of the first BZ (0, 0, 0.5) during the frequency range from about 2 to 3 THz (only cyan bubbles there).

Conclusions and discussion
To conclude, the thermal conductivity in silicon nanowires varies from the center to the surface. It is maximum at the center and decreases close to the surface. It is observed with both MD and MC simulations. For diameters smaller than the dominant intrinsic phonon MFP, the radial evolution of thermal transport does not depend on the total diameter of the nanowire. At macroscopic scale, the reduction of the thermal conductivity close to the surface vanishes and is not visible anymore. The reduction of the heat flux at the surface can be due to surface scattering and/or modifications of phonon properties close to the surface (amorphous-like local DOS, phonon localization, K). The appearance of localized surface modes could also be at the origin of the decrease of the thermal conductivity. The combination of these different phenomena could lead to even further reduced thermal transport close to the free surfaces in all semi-conductors nanostructures.
During experimental measurements of the thermal conductivity in nanowires, a special attention should be given to the radial evolution of heat transport. With the 'bridge' method, where the nanowire is suspended between a hot and a cold bath, the measured axial thermal conductivity is the total axial one. However, when using a tip to heat the surface of a nanowire, the measured heat current depends first on the radial conductivity and then on a average of the axial and radial one. If the nanowire is only heated close to the surface, the measured Figure 6. Phonon bands (light gray solid line) of the Si NW with d=5.4 nm which is composed of 5 concentric cylindrical layers with the same thickness 0.54nm. The geometric structure with layers from internal to external sequentially labeled by the numbers from 1 to 5 is shown on the right. For different layers, the average partial contribution of each atom in each phonon mode are overlaid on the total phonon dispersion. The bubble sizes of the bands are scaled by the amplitudes of specified vibrations. For the sake of clarity, only modes with average partial contribution larger than 0.35% are overlaid by bubbles in different colors corresponding to different layers as indicated on the right. conductivity is not the axial one. This could explain the discrepancies between thermal conductivities given by different experimental methods observed in nanofilms or nanowires. When measuring the thermal conductivity of nanowires with the tip method, if the penetration depth is small compared to the diameter, a correction factor of at least 1.2-1.4 (MC-MD) should be applied to get the total thermal conductivity along the nanowire axis.