Effect of Boundary Reflectivity on Thermal Transport Properties of Single-Walled Carbon Nanotubes Examined by Molecular Dynamics Simulations

We study the effect of tube-end reflectivity on thermal transport property calculations of carbon nanotubes with non-equilibrium molecular dynamics to explore the possibility to reduce the length of heat bath which has to be comparable to the length of the region of interest in the case of rigid boundaries (RB) adopted in previous studies. We implement atomic scale non-reflecting boundaries (NRB) for this purpose, and found that the NRB acts as a thermal sink, which is not ideally suited to the task of improving thermal property calculations; however, a hybrid asymmetric system with the NRB condition for one tube-end and the RB for the other is found to reduce cold end resistance to near zero, and thus is promising to improve the efficiency of thermal property simulations while keeping the accuracy. [DOI: 10.1380/ejssnt.2010.313]


I. INTRODUCTION
The continued decrease in the size of electronic devices has lead to the need for improved heat transfer in nanoscale devices.This is due in part to the increased current density required to transmit data through the reduced component size, which results in increased joule heating and leads to defects such as thermo diffusion [1,2].Carbon systems have been shown to exhibit unique and interesting thermal properties [3], which gives hope to further component size reduction for greater speed.
The high thermal conductivity of single walled carbon nanotubes (SWNT), comparable to that of diamond (about 2000 W/m K) [4][5][6], have made them ideal candidates for thermal transport in nanoscale technology in recent years.It's interesting to note that the chirality, radius, and length of these quasi-one-dimensional conductors can change thermal transport properties of SWNTs [3,7].Furthermore, carbon nanostructures can be further categorized into a variety of novel geometric features with unique thermal properties.For example, thermal rectifying properties have recently been predicted theoretically in carbon nanocones [8] and asymmetric graphene ribbons [9], which may lead to realizing future phononics applications.It has also been shown that randomly arranged SWNTs can curtail the high conductivity and act like an insulator [10].These nanostructures are promising candidates for future thermal transport applications and give hope to achieving phononics applications in the near future.
In spite of the promising outlook, the response of experimental systems to defects (e. g. vacancies, adatoms, probe, etc.) is not well understood and may have undesirable consequences on more complex, large scale systems.Simulations are powerful tools, and can be used to study the effect of small scale defects on large scale thermal transport systems.However, carbon nanostructures need integration into a multiscale thermal transport environment in simulations for this purpose.
Current atomic-scale simulation techniques have an- other subtle problem, that is, thermal boundary resistance (TBR) between reservoirs and the freely conducting region.Figure 1 shows a typical atomic configuration of molecular dynamics (MD) simulations for thermal properties and schematic spatial temperature profile obtained from simulations using the configuration adopted here and in previous simulations of this kind.In the configuration, the region of interest, where temperature is uncontrolled, is sandwiched between cold and hot reservoirs, and rigid boundary (RB) conditions are imposed at the tube-ends.As shown in the temperature profile of Fig. 1, temperature drops corresponding to TBR appear at the interfaces between the reservoir and the uncontrolled region when the length of temperature-uncontrolled region is less than phonon mean free path.It is known that the non-equilibrium MD simulation under the RB condition overestimates the temperature drop when the length of each reservoir, L R , is much shorter than that of the temperature-uncontrolled region, L B .This originates from unphysical phonon reflection at the boundary.Previous studies have shown that the length of each reservoir is needed to be the same as or more than the length of the uncontrolled region to remove the artificial phonon reflection [11].This poses two limitations on continuing to multiscale thermal transport simulations: 1) computa-tional inefficiency due to large reservoir size and 2) nonideal, resistive transition between reservoir and conducting regions in the atomic-scale space.It would be ideal to simulate diffusive transport using atomic scale systems of size generally dominated by ballistic or quasi-diffusive transport.
This study investigates the effect of reflectivity at tube ends as a way of controlling thermal boundary resistance.We implemented a molecular non-reflecting boundary (NRB) condition [12,13], which is the first ever use of such a boundary on carbon systems to our knowledge.The NRB method proposed earlier has been used previously as a pressure transmitting boundary to reduce shockwave reflection in crystalline systems [14,15].Here we apply the first ever application to atomically thin systems.The lack of periodicity and crystallinity in the carbon nanotube allows us to investigate application of further degrees of freedom to this NRB in future studies such as radial and torsional NRB; however, we will focus on longitudinal NRB in this study.We present here the first steps to producing a NRB condition for carbon molecular systems.Both elastically reflective rigid boundary (RB) and longitudinally NRB are investigated for their effect on thermal transport systems.In addition, results of an asymmetric boundary, NRB for one tube-end and RB for the other, shows promise as a way of reducing computational size along with TBR.

A. Thermal transport simulations
Phonon dominated heat transport has been observed in SWNTs at room temperature theoretically [16,17] and experimentally [18], making non-equilibrium molecular dynamics (MD) simulations the ideal tool for this study.Molecular dynamics is therefore adopted by this study along with the Brenner interatomic potential [19] because of its reproduction of thermal, elastic, and structural properties of carbon.
The Nose-Hoover thermostat is used because it ensures a canonical distribution [20,21].For SWNTs, a temperature difference, ∆T , of 20 K is used between reservoirs held at 290 K and 310 K for left and right end temperatures.
Thermal boundary resistance is quantified as R L and R R , the resistances of the left and right reservoirs respectively such that R i = (linear fitted temperature drop) / ∆T , where the linear fitted temperature drop is defined as the difference between the reservoir temperature and the temperature at the boundary estimated a by linear regression of the temperature-uncontrolled region.Note that the ideal temperature gradient, as measured longitudinally along z, is dT I /dz = ∆T /L B such that TBR is zero, corresponding to completely diffusive transport independent of L B .Realization of the ideal case would allow for simulations of longer nanotubes connected to the thermostat by continuum or mesoscopic models with resistance that is more physical.
The non-dimensional parameter η = L R /L B is used as a measure of relative reservoir length to bulk length (see Fig. 1(b)).A (5,5) nanotube with length L B = 25 nm is used in this study, which is long enough to show quasi-diffusive transport but very nearly ballistic.It is the hope of this study to use this short nanotube to show the ability to achieve improved results in nearly ballistic systems.Note that dT I /dz= 0.8 K/nm when L B = 25 nm and ∆T = 20 K.
The physical quantities shown are all averaged over 2 ns with a time step of 0.5 fs.The time interval used is sufficient for all simulations to result in steady state thermal current (J th ) conditions to be satisfied.J th is calculated in the form of Kondo et al. [22] assuming a cross-sectional annular radius equal to the time-averaged radius of each case (∼3.49Å) and thickness 3.4 Å, the interlayer spacing in graphite.Before the temperature gradient was applied, all SWNTs were brought to equilibrium at 300 K for 500 ps.

B. Non-reflecting boundary (NRB) condition
It is thought that one cause of TBR in simulations is the unphysical reflection at the tube end.Every atomic layer in the nanotube has some finite probability to reflect an incoming wave.The probabilities of reflection off of a RB and NRB are, in principle, 1 and 0 respectively.In principle the probability is less than one for a nanotube that extends beyond the simulated boundary.
At the tube-end boundary, the force acting on an atom is different from the one in the realistic (ideal) system without boundaries.The ideal force on a boundary atom F B (i.e. the force expected in the ideal system without boundary) can be calculated as Here F MD is the force in the z-direction obtained from MD simulation and F NRB is the force from the atoms on the other end of the boundary which are removed from the system to make the number of atoms finite.As an example, consider any pristine atomic layer of the carbon nanotube in Fig. 1.In the case that we wish to reduce the number of simulated atoms by removing all atoms to the left of our selected atomic layer, we then define F B such that F MD is the force of all MD simulated atoms to the right of our layer and F NRB is the estimate of the removed atoms.The opposite set of forces would be used if atoms were removed to the right of the atomic layer.It should also be noted that because the atomic layer is pristine, no interlayer force is in the z-direction; therefore, we neglect all interlayer forces.In the NRB condition, the forces on the boundary atoms are set to F B in eq.(1).To estimate F NRB , the following linear form had been proposed [12,13]: ( Here, F 0 represents the z-direction force required to keep the boundary at its equilibrium position and ζ is the characteristic impedance multiplied by cross-sectional area for an atomic layer with instantaneous longitudinal velocity V z per atom.Under these conditions, the NRB as proposed earlier and used here can be viewed as a onedimensional damper that reacts to remove longitudinal waves impacting the boundary.eq. ( 2), we examined the relation between V z and F MD from atoms of lower z (i.e.located on the left side) for atoms inside the simulation region except the tube ends.For these atoms, F MD from atoms of lower z, which corresponds to the exact F NRB in eq. ( 1), is available from MD simulations.The result is shown in Fig. 2, and as we can see from the figure, clear linear relationship is seen.Though the values are scattered and thus correlation is not perfect, we can say that the use of eq. ( 2) is justified as an approximation.In practice, F NRB was estimated using a 20 nm (5,5) SWNT with a 3 nm longitudinal square velocity wave applied at 300 K.After several simulations, F 0 and ζ were averaged to approximate 0.84 eV/ Å per atom and 0.123 eV-ps/ Å2 respectively.On the other hand, in principle, ζ can be estimated by the product of density ρ, speed of sound c, and cross-sectional area S to be approximately 0.19 eV-ps/ Å2 per atom in the atomic layer.
Using the square impulse wave, c was found to approximately equal 19 km/s.The agreement on the value of ζ between this estimate and the one by MD simulations is satisfactory.
As seen at the tube ends in Fig. 1, two atomic layers represent a boundary.For the NRB condition, each atomic layer is an independent NRB such that F MD is the force of atoms towards the center of the tube and F NRB is the linear estimated from eq. ( 2).It should be noted that the force of an outer layer is not directly applied to the inner layer.Because the Brenner potential is significant to the second neighbor shell, it becomes important to have the motion of the outer layer for the atoms adjacent to the boundary.
In addition to the system with the NRB condition, an asymmetric system having the NRB condition at the left tube-end (cold) and the RB at the right end (hot) was investigated.This was done in order to see effectiveness of removing the coherency of the vibrating tube ends.Furthermore, asymmetry of the boundary reflectivity will impose preferential conduction towards the NRB, situated at the cold end in this study.
In a realistic carbon nanotube, phonons created both inside the MD region and in the nanotube beyond the boundaries can enter and return.However, the NRB lacks the ability to bring waves back into the system from the background SWNT.It is hypothesized that the asymmetric boundary may represent return of phonons from beyond the RB into the system in some fashion.We therefore use it here as a first step to understanding the effect of asymmetric boundaries, though the real experimental situation corresponding to such boundary conditions is unclear.

III. RESULTS AND DISCUSSIONS
The three systems NRB, RB, and NRB-RB were each tested for their effect on the time-averaged temperature gradient within the uncontrolled region for a variety of η, η= 0.05, 0.15, 0.25, 0.5 and 1.In order to reduce statistical fluctuations, gradient and resistance properties are shown as the average of two simulation runs.As can be seen in the obtained temperature gradients (Fig. 3), the systems with a NRB have reduced mean temperatures compared to RB with little effect on the temperature gradient.This is as expected because the boundary annihilates most incoming longitudinal waves, which acts as a form of thermal sink.As expected, the temperature gradient increases with η for each boundary case except minor reduction in the RB case for η = 0.15 to 0.5, which is clearly seen in Fig. 4(d).
Figure 4 shows the thermal properties, R L , R R and thermal current as well as temperature gradient, of the three boundary cases.The total cooling in the NRB case results in negative resistance for R L in the NRB and NRB-RB conditions.For NRB, the total contribution of resistance to the temperature gap surpasses 70% of reservoir difference in the η = 0.05 simulation.The difference in total TBR (|R L | + |R R |) from the RB case is statistically negligible for η > 0.25 (Figs.4(a), (b), and (c)).However, it is important to note that R L fluctuates around 0 for η > 0.25.This may be the result of cooling in the entire system resulting in a total mean temperature reduction equal to the left end temperature drop.This result is promishttp://www.sssj.org/ejssnt(J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) ing when seen combined with the thermal gradient (Fig. 4(d)) of the NRB-RB system, which is greater than both RB and NRB if compared at the same η.As mentioned in the introduction, saturation of thermal gradient is expected for sufficiently large η.In fact, our results in Fig. 4(d) exhibit such a tendency for all the three boundary condition cases.The behavior of the NRB-RB mentioned above suggests that thermal current reaches the saturated value with the increase in η in this case faster than the RB and NRB systems.This means that we can expect to achieve a more ideal thermal gradient simultaneously with one end of the thermal gradient having near seamless thermal transition to the thermostat.The asymmetric case is therefore promising to achieve our goal of zero TBR and ideal gradient.It was also observed that J th was significantly reduced in the NRB-RB case (see Fig. 4(e)).It was found that when both reservoirs are held at 300 K, a thermal current appears.For symmetric boundaries, the same conditions results in J th = 0.The cause of the background current is likely the result of local cooling at the cold reservoir, which causes an increase in energy to the cold thermostat compared to the hot reservoir.Therefore, for non-symmetric boundaries, we a background cur-rent (J B ), which is of the opposite sign of the boundary induced current (see Fig. 5), which increases with the size the system.Therefore, the thermal current due to the temperature difference of an asymmetric SWNT is J = J th − J B .Applying this correction, we get the teal curve in Fig. 4(e), which is comparable to the symmetric NRB and RB conditions.
The presence of background current in asymmetric boundaries shows that there may be a small boundary induced temperature gradient between asymmetric tubeends.The NRB case may appear to represent a more realistic system because it equally transports incoming waves out of both ends of the nanotube.However, a true conducting nanotube will have phonons not only leaving, as in the case of NRB, but also returning or being formed outside of our narrow MD region and entering through the simulated tube end.Additionally, in the event that a temperature drop occurs between both reservoirs, there would be a reduction in total energy of waves entering through the cold end relative to the hot end.In this case, the choice of boundary conditions affects the interpretation of conduction in the system beyond the boundaries.The true nature of the asymmetric boundary cannot be fully justified here nor is the effect on the phonon spectra the subject of this study.However, the observed improvements to the thermal transport simulations make asymmetric boundaries a direction to focus future investigations.

IV. CONCLUSIONS
Non-equilibrium MD simulations of SWNTs were performed with boundary conditions of varying reflectivity.Non-reflecting boundaries act as a thermal sink, which lowers the overall temperature of the uncontrolled region.This has the unintended consequence of reducing the resistance at the cold reservoir for larger reservoir sizes, which is ideal.This improved resistance along with observed increases in thermal gradient for the asymmetric boundary condition make it the best case for further development of these techniques because it acts as a system larger than the mean free path while remaining computationally efficient.
The asymmetric case also was found to have a background current induced by the presence of asymmetry in the boundary condition.The exact nature of this consequence on the system and how it improves thermal properties is not well understood and will be the focus of future research.Further development of this method, including investigation of symmetric and asymmetric partially reflecting boundaries, will certainly be an area of future research.It is curious that all of these cases caused reduction of resistance at the cold boundary, but further investigation is needed to reduce resistance at the hot reservoir.
As the resistance gap separating the atomic scale mod-eling technique to the continuum thermostat is improved, these systems will be more efficient and capable of handling smooth transition between models.This is just the next step towards reducing the computational error in these techniques to make future theoretical models more precise.

FIG. 1 :
FIG. 1: Schematic representation of heat transfer simulation showing the spatial temperature profile with thermal boundary resistance (TBR) and the atomic configuration.Tube ends (grey), cold reservoir (blue), hot reservoir (red), and uncontrolled region (green) are shown.

FIG. 3 :
FIG. 3: Temperature gradients for three boundary condition systems.Black arrows show the effect of increasing η on the fitted line for RB (a), NRB (b), and NRB-RB (c) conditions.