Effects of ballistic transport on the thermal resistance and temperature profile in nanowires

Abstract Effects of ballistic transport on the temperature profiles and thermal resistance in nanowires are studied. Computer simulations of nanowires between a heat source and a heat sink have shown that in the middle of such wires the temperature gradient is reduced compared to Fourier’s law with steep gradients close to the heat source and sink. In this work, results from molecular dynamics and phonon Monte Carlo simulations of the heat transport in nanowires are compared to a radiator model which predicts a reduced gradient with discrete jumps at the wire ends. The comparison shows that for wires longer than the typical mean free path of phonons the radiator model is able to account for ballistic transport effects. The steep gradients at the wire ends are then continuous manifestations of the discrete jumps in the model. Graphical abstract


Introduction
The development of nanotechnology and the miniaturization of devices has lead to a strong interest in the subject of thermal transport at the nanoscale.As is the case with other physical properties, the laws describing thermal transport in macroscopic systems must be modified on the nanoscale.This is not only important for the design of nano-devices but can also be exploited to create nanomaterials with unusual thermal transport properties.Examples of technological areas that can benefit from an understanding of nanoscale thermal transport can be found in [1,2].
Heat transport in macroscopic solid systems is well described by Fourier's law which is based on the idea of a diffusive transport mechanism where heat carriers are Brownian particles which undergo frequent collisions that lead to random direction changes.However, when the system dimensions are comparable to or shorter than the mean free path of the heat carriers between collisions, Fourier's law no longer provides an adequate description and microscopic transport mechanisms must be considered.
In this work, we are using the simple case of a nanowire between a heat source and a heat sink to study deviations from Fourier's law caused by ballistic transport effects.Results obtained from molecular dynamics simulations and phonon Monte Carlo simulations are compared to the prediction of Fourier's law and a simple phonon radiator model.
Figure 1 shows the type of systems we are considering in this work.A nanowire is connected to two thermal heat baths with temperatures T l and T h .Fourier's law predicts that between the heat cold lead hot lead nanowire baths, the temperature follows a straight line that connects the temperatures of the two heat baths.Several studies of thermal transport through nanowires or thin films have shown deviations from the behaviour predicted by Fourier's law.Solutions of phonon radiator models [3,4] predict a linear temperature profile with a reduced slope and discrete temperature jumps at the interface to the heat baths.A similar behaviour can be seen in results from a phonon Monte Carlo simulation (Fig. 4 in Ref. [5]).Molecular dynamics simulations [6][7][8][9][10][11] on the other hand show non-linear temperature profiles with steep slopes near the areas that are thermalized through contact with a heat bath (leads).Away from the leads the temperatures in the molecular dynamics simulations settle into a linear behaviour with a reduced slope compared to the Fourier case.
The deviations from Fourier's law observed in the previous studies can be understood qualitatively by looking at the insets in Fig. 1.The temperature profile in the immediate vicinity of any point in the middle of the wire (right inset) is symmetric so that the temperature increase to the right of the point is compensated by a similar decrease to its left.Further, the shape of the nearby profile remains the same everywhere in the middle of the wire.The situation changes, however, when we move the point close to one of the leads.As shown by the right inset, the temperature profile in the immediate vicinity is no longer symmetric and its shape depends on the distance to the lead.Since the temperature profile of the surroundings determines the distribution of incoming heat carriers, the behaviour of the temperature close to the leads has to account for the change in the environment.
The goal of this work is to show that the reduction of the temperature gradient inside the wires is a general phenomenon that is well described by a simple radiator model.To this end, we compare the analytical predictions of the model to results from molecular dynamics and phonon Monte Carlo simulations.This comparison shows that the radiator model provides a quantitative description of the reduction of the temperature gradient over a wide range of wire lengths.
The rest of the article is organized as follows.In Sec. 2 we present the theory of the radiator model.(In the appendix we also present an extension of the model to a periodic geometry which facilitates the comparison to the molecular dynamics simulations).In Sec. 3 and Sec 4 we compare the predictions of this model with results from molecular dynamics and phonon Monte Carlo simulations, respectively.The article finishes with a brief summary and conclusion in Sec. 5.

Model calculations
To understand the behavior of the temperature close to a thermalized lead, we consider a simple model for heat transport based on the idea introduced by P. Debye in Ref. [12].In this model, some kind of heat carriers store and transport heat energy by flowing through the material while undergoing random scattering events.The model makes no assumptions about the nature of the heat carriers and does not take into account different types of carriers such as phonons with different wavelengths or polarizations.A discrete version of a similar model was used by Klitsner et al. [3].
For a thin nanowire, we can treat the system as a quasi-one-dimensional problem so that there are only two possible directions for heat or energy flow along the wire axis.For simplicity, we refer to these two directions as left and right, independent of the orientation of the wire.For a thin slice of the wire at position x we thus have to consider only the two energy fluxes Φ l (x) and Φ r (x) entering the slice from the left and right side, respectively.
Scattering reduces the energy fluxes as they pass through the slice.Assuming a Lambert absorption law the flux reductions are proportional to the magnitude of the fluxes and the width of the slice ∆x.With an attenuation coefficient µ, the total rate of energy going into scattering events, S, in the slice at position x is then In the steady state there can be no net gain or loss of energy unless the slice is attached to a heat source or sink.Energy scattered from the incoming fluxes by scattering events is therefore re-emitted from the slice by adding the amount 1 2 S(x) to the fluxes leaving the slice on either side.
If the attenuation coefficient µ is the same everywhere in the system, the Lambert absorption law leads to an exponential decrease of energy fluxes travelling through the system.Concretely, if 1  2 S(x ′ ) is the flux emitted by a slice at x ′ in one direction, only the amount reaches a slice located at x. Equation 2 allows us to calculate the total fluxes Φ l (x) and Φ r (x) entering the slice at x.These two fluxes can be seen as the sum of contributions of an infinite series of slices to the left and right of the the slice, respectively.Taking the continuum limit of infinitesimally thin slices, and replacing the summation by integration, we obtain: Combining these expressions with Eq. 1 we obtain the steady state equation To calculate temperature profiles, we need to connect the rate of energy being scattered by the heat carriers S with the temperature T .Debye argued that the fluxes passing through a crosssection of the wire is proportional to the internal energy density of the material [12].In the high temperature limit -well above the Debye temperature -the internal energy is proportional to the temperature T .Since the rate of energy being scattered at x is proportional to the fluxes at this point, we assume that S(x) = αT (x).Putting this into Eq. 4 and removing the proportionality constant α that appears on both sides, yields This steady state equation determines the temperature at locations with no heat source or sink.
For a wire of length L, connected to infinitely long thermal leads with constant temperatures T l and T h , it can be shown that Eq. 5 is solved by a linear profile A proof of this solution is given in the appendix.Our simple model thus predicts a linear temperature profile similar to Fourier's law.However, the temperature gradient predicted by Eq. 7 is lower than in case of Fourier's law where the gradient is T h −T l L .The lower gradient results in discontinuous temperature jumps at the interface between the leads and the wire.Figure 2 shows the temperature profile for various values of the attenuation constant µ.The size of the temperature jumps decreases as the length of the wire L increases.In the limit of L → ∞ the jumps vanish and Eqs.6, 7 converge towards the prediction of Fourier's law.
The discontinuous temperature jumps predicted by the model are the necessary accommodations described in the introduction.These jumps decrease the fluxes generated by the linear profile inside the wire so that they equal the fluxes entering from the constant temperature leads.Without  these jumps, the change from the flat temperature profile in the thermalized areas to the linear profile in the middle part of the wire would lead to an inconsistency between the temperature fluxes entering the system from the leads and the flux inside the wire.We continue the discussion of the discrete temperature jumps at the end of this section.
Further insight into the effect of the leads can be obtained by calculating of the net flux inside the wire as the difference between the fluxes going in the left and right direction: Using Eqs.6,7 we obtain for −L/2 < x < +L/2 The net flux is thus constant throughout the wire as it should be since there is no heat source or sink in this area.From Eq. 9 we find the thermal conductance of the wire system The thermal resistance is easier to interpret.Equation 11shows that the resistance consists of two contributions.The first term is the thermal resistance of the wire material.It is proportional to the length L of the wire with a bulk conductivity of κ ∞ = µ α .The second term represents a constant resistance R c = 1 α that is added at the interface between the wire and the leads.Note that this is not the usual Kapitza resistance that occurs at interfaces between different materials.Our model assumes the same material for the lead and the wire.The only difference between the wire and the leads is that the leads are kept at constant temperatures.
In summary, the simple model presented in this section predicts the occurrence of contact resistances at the interfaces between the constant temperature leads and the linear temperature profile inside the wire.These contact resistances cause discrete jumps in the temperature at the interfaces.They reduce the net flux inside the wire so that it is consistent with the fluxes injected by the leads.This model is of course a simplified picture that neglects many details of the transport process such as different phonon wavelengths and polarizations.We expect that a more realistic treatment of the transport process would result in a smooth temperature profile without discrete jumps.

Molecular dynamics simulations 3.1 Details of the simulations
In this section we present results from a series of non-equilibrium molecular dynamics simulations of nanowires.The goal of these simulations was to understand to what extent the model calculations presented in the previous section predict the outcome of a more realistic treatment of the thermal transport on the nanoscale.Since molecular dynamics simulations are limited to finite system sizes, the simulations were carried out with finite simulation cells using periodic boundary conditions along the wire axis x.Along the y and z directions open boundary conditions were applied.As shown in Fig. 3 mirror symmetry of the simulation cell, the temperature profile in the right half has to be a mirror image of the profile in the left half.In our analysis we therefore averaged the profiles of the two halves.As pointed out in the appendix, the model discussed in Sec. 2 can be also be solved for the periodic geometry used in the molecular dynamics simulations.The slope of the profile in the wire section for this geometry is given by Eq.B8.
The simulated wires were built by cutting a cylinder with a radius of 10 lattice constants from a perfect fcc lattice.The cylinder axis was aligned with the lattice's [100] direction.For the interaction potential we chose the Lennard-Jones potential [13].To keep our results independent from any specific material, we set the energy and length parameters ϵ and σ of the potential as well as the particle mass m to unity.All parameters and results of the simulation are presented here using these dimensionless (reduced) units.All simulations were carried out with the simulation package LAMMPS [14] using a time step ∆t = 0.005.
Our first step was the determination of the equilibrium lattice constant of the wire systems along the wire axis.This step was necessary since the finite temperature as well as relaxations of the crystal lattice perpendicular to the wire axis lead to deviations of the lattice constant from its value at T = 0.The lattice constant was determined from a simulation of a wire with a length of 200 lattice constants.During this simulation, the size of the system was relaxed with the help of the Parrinello-Rahman method [15,16].The lattice constant a 0 = 1.55030 determined from this simulation was then used in the following non-equilibrium simulations.
The main part of our molecular dynamics simulations was a series of simulations using wire systems of different length.The wire lengths in these simulations were: L = 25 a 0 , 50 a 0 , 100 a 0 , 200 a 0 , 300 a 0 , 400 a 0 , 500 a 0 , 600 a 0 , 700 a 0 , 800 a 0 , 900 a 0 , and 1000 a 0 .The width of the thermal leads was in all instances w = 50a 0 .The total number of atoms in these configurations ranges from 187,350 to 2,622,900 atoms.The thermal leads were realized in these simulations with the help of two Langevin thermostats set to temperatures T l = 0.045 and T h = 0.055.The damping constant of the thermostats was set to 1 time unit.
For the computation of temperature profiles in the simulated wires, we subdivided the configuration into slices with a thickness of 2.5 a 0 along the wire axis.During the simulations, the number of atoms in each slice and their total kinetic energy were measured each tenth simulation step and accumulated over periods of 1000 steps.At the end of each 1000 step period, the equipartition theorem was used to calculate a local temperature for each slice during the period.These data were written into a file which was used to monitor the progress of the simulation and to calculate profiles averaged over longer time spans.
Prior to the measurement of the temperature profile, we simulated all systems until they had reached their steady state.In addition to a manual inspection of the temperature profiles, we used the time-evolution of the average slope of the wire section as the main criterion to decide when a system had reached the steady state.To this end, we averaged the temperature data over periods of 100,000 simulation steps and applied a linear fit to the temperature data in the wire section.The simulations were continued until the slope of the fit line no longer showed a systematic change and had settled into fluctuations around a constant value.

Temperature Profiles
After the model configurations had reached the steady state, the simulations were continued for another 1,000,000 steps.From the temperature data written to file during these simulations, the average temperature profiles of all configurations were determined.As an example, we show in Fig. 4 the profile calculated for the system with L = 100 a 0 .The figure shows a very good agreement between the simulations and the predictions of the model presented in Sec. 2. There are sharp temperature jumps at the interface between the thermal leads and the wire.Apart from a small transition region close to these interfaces, the profile follows a straight line inside the wire and remains constant in the thermal lead sections.The temperature gradient in the simulation (green line) is clearly reduced compared to the prediction of Fourier's law (thin black line).A similar behaviour was observed for the systems with other wire lengths L.
The presence of a transition area close to the interface between the thermal leads and the wire is to be expected for the same arguments we gave in the introduction.The thermal flux in the wire constantly adds or removes energy from the outer parts of the thermalized leads so that the environment there is different than in the centre of the leads.The thermostats can only heat or cool the leads at a certain rate, which depends on the thermostat settings.The temperature profile near the border of the leads is the result of a balance between the flux delivered by the thermostat and the flux going into or coming from the wire.A similar result might be expected in an experimental situation where the contact to a heat bath provides only a limited heating or cooling rate.
In order to study the effect of the thermostat settings on the temperature profile, we have repeated the simulations of the wire of length L = 100 a 0 with the Langevin damping parameter scaled up and down by a factor of 4 as well as a Nose-Hoover NVT thermostat.The results of these simulations are given in Fig. 5.In order to show the differences between the setting clearly, we only show the transition at the lower temperature end of the system in this figure.The results at the higher temperature end are qualitatively the same.
The temperature profiles for different Langevin damping parameters in Fig. 5 show that lower damping parameters result in a sharper temperature gradient at the interface.This agrees with the mechanism described in the last paragraph since lower damping parameters allow for a faster energy transfer between the thermostats and the systems.
Figure 5 shows that the NVT thermostat behaves different than the Langevin thermostats.Apart from a very smooth entry into the thermostatted area, the temperature undershoots the set temperature of the thermostat (0.045) in the middle of the thermostatted region.This can be understood from the mechanics of this thermostat.The NVT thermostat does not drive the temperature (more precisely its kinetic energy) of each atom towards the set value.Instead, the thermostat maintains the average temperature in the thermostatted region.The higher temperatures close to the interface in Fig. 5 are therefore balanced by lower temperatures in the centre of the thermal lead.
While Fig. 5 shows that the thermostat parameters affect the profile in the thermal leads, the figure also shows that there is little variation between the profiles inside the wire.Most importantly, the temperature gradient in the middle part of the wire is virtually the same in all four cases.This indicates that the deviations from Fourier's law are not caused by the dynamics of the thermostats.
In order to make a quantitative comparison between our simulation results and the model predictions, we determined the slope of all profiles with the help of a linear regression of the temperature data in the wire sections.In order to minimize the effect of the transition region at the interfaces to the leads, the first and last data points were excluded from the fit.In the left panel of Fig. 6 we present the resulting slopes or temperature gradients as a function of the wire length with the predictions of Fourier's law.It can be seen that the temperature gradients in the simulations are systematically below the prediction of Fourier's law (orange line) but converge towards the latter for longer wires.This is the expected behavior since Fourier's law is correct for macroscopic systems.Further details are revealed by the inverse slope data shown in the right panel of Fig. 6.It can be seen that for longer wires the inverse slopes follow a straight line parallel to Fourier's law.This is in agreement with Eq.B8 which predicts the inverse slope to be The blue lines shown in Fig. 6 were obtained by fitting the constant term in Eq. 12 to the simulation data for L ≥ 200 a 0 .It corresponds to an attenuation constant µ = 0.0684.It can be seen from Fig. 6 that Eq.B8 (represented by the blue fit line) provides a better description of the temperature profiles resulting from the simulations than Fourier's law.For shorter wires, however, Eq.B8 underestimates the temperature gradients.In our opinion, this discrepancy is related to the fact that the temperature profiles in the simulations do not show a discrete temperature jump as predicted by the model in Sec. 2. Instead, they feature a thin but finite width transition region between the thermal leads and the wire (see for example Fig. 4).Since the transition regions extend into the thermal leads, they reduce the need for adaptation in the wire section.Such a transfer of the adaptation is excluded in the model calculations which keep the lead temperatures constant as a boundary condition.The impact of the finite transition areas would be even larger for shorter wires since the temperature jump increases for shorter wires.
Another factor that might contribute to the differences between Eq.B8 and the simulation data is our method to compute the slope values from the simulations.We carried out a number of tests where we changed the number of points excluded from the linear fit of the wire section.We also experimented with some non-linear fit functions that could account for the transition regions.None of these attempts lead to significant changes of our results.We are therefore convinced that the presence of the transition areas which is not included in the model calculations is the main reason for the increased temperature gradients in the simulations of shorter wires.
Overall, we believe that the results presented in this subsection confirm our basic premise that the temperature profile of the wire reacts to the different environment near the thermal leads.In accordance with our model calculations, the resulting accommodation at the wire ends result in a notable reduction of the temperature gradient at the center of the wire.

Thermal resistances
In addition to the temperature profiles, we used the energy flux in the wires as a second quantity to compare the predictions of our model calculations to the molecular dynamics simulations.During the simulations, the cumulative amount of energy added to or subtracted from the system by the two thermostats that maintain the temperatures T h and T l in the lead sections were computed and periodically outputted (together with the temperature data).The difference of the two tallies corresponds to the total amount of energy that has been transported through the wire sections during the simulation.Over the course of the simulation, this value increases on average linearly with time.By fitting a straight line to the data and taking its slope, we were thus able to determine the average net energy flux running through the wires and from these the thermal resistance R for each wire configuration.
In Fig. 7 we present the thermal resistances R determined from the molecular dynamics simulations as a function of the wire length L. This figure  agrees qualitatively very well with Eq. 11 which predicts that the thermal resistance is a linear function of the wire length with a non-zero offset.A complete quantitative comparison between the predictions of Eq. 11 and the simulation data is not possible since the model does not yield information about α.However, according to Eq. 11 the attenuation constant µ equals twice the slope of the fit line divided by its offset.Applying this method to the values of the fit line in Fig. 7 gives a value of 0.0532a −1 0 which is about 22 % smaller than the value µ = 0.0684a −1 0 obtained from the temperature profiles.
In our opinion, the difference between the µ values is a result of the fact that the model used to derive Eq. 11 does not account for the details of the transport mechanism.We believe, however, that Fig. 7 provides support for the general prediction of the model which is that the interfaces between the thermal leads and the wire give rise to a constant interface resistance which induces steep temperature changes in the vicinity of the interface.It should be noted that the results in this subsection are independent from the results obtained in the preceding subsection since the cumulative energy tallies of the thermostats are not affected by the local temperature distribution inside the wire.

Phonon Monte Carlo simulations 4.1 Details of the simulations
The results presented in the preceding section show a good agreement between the molecular dynamics simulations and the predictions of the model from Sec. 2. In this section we present results from phonon Monte Carlo simulations as an independent test of the validity of the model predictions.By using a completely independent simulation method we are able to remove any remaining doubt whether the results in the preceding section are artefacts of the artificial dynamics of the thermostats in the molecular dynamics simulations.
The phonon Monte Carlo simulation method is a relatively new method for the simulation of heat transport in non-metallic materials.The method was first introduced by Mazumder and Majumdar [5] based on preliminary work by Peterson [17].Further improvements were published by Lacroix, Joulain and Lemmonier [18] as well as Péraud, Hadjiconstantinou and Nicolas [19,20].In a phonon Monte Carlo simulation, packages of phonons are propagated through the simulated system.During the simulation, the phonons in the packages undergo stochastic collisions leading to changes of the direction and momentum (frequency) of the phonons.Details of the method can be found in Refs.[5,[18][19][20][21].Details of the implementation of the phonon Monte Carlo method used in this work can be found in [22].
Phonon Monte Carlo simulations require a scattering model that describes the lifetime of different phonon modes.In this work we have used two different such models: the more recent parametrization used by Jean et al. [23] (Si-1 ) and the older parametrization used by Lacroix et al. [18] which is based on data by Holland (Si-2 ).The usage of the two models allows us to study the effect of the phonon mean free path since the latter model is known to predict shorter phonon life times and therefore shorter mean free paths than the former.Both models only account for the phonon scattering rates due to phonon collisions.It would be possible to add other scattering mechanisms such as impurity scattering and/or partially diffusive surface scattering to the simulations.These effects, in particular diffusive surface scattering, are certainly important for a quantitative description of thermal transport in nanowires.However, for the purpose of this work which is to study the effects of the thermal leads, the difference between Si-1 and Si-2 is sufficient.A more detailed study of the effects of different scattering mechanisms is beyond the scope of this work.
Using the two scattering models, we performed simulations of straight wires with lengths L in the range from 1 µm to 50 µm in the case of Si-2 and 0.1 µm to 10 µm in the case of Si-2.The wire width was in all cases 20 nm.It should be noted that the simulation code only allows the specification of a two-dimensional geometry.The z-dimension is treated as an infinite open system.Apart from a scaling factor, the results obtained in this manner are, however, equivalent to a system that is limited along the z direction by perfectly specular boundaries.Our results therefore correspond to simulations of a square wire with 20 nm side length.For all wire lengths, we simulated 500,000,000 phonon packages.with a maximum lifetime of 400 ns per µm of system length, L.

Temperature Profiles
As an example, Fig. 8 shows the temperature profiles of a Si wire with a length L of 2 µm obtained from phonon Monte Carlo simulations for the two scattering models.It should be noted that the phonon Monte Carlo simulations do not simulate the thermalized leads but consider only phonons emitted from these areas into the wire (or absorbed by them).The lead areas in Fig. 8 are therefore not simulation results but are shown in order to visualize the change of the profile, most notably the jumps at the interfaces.
The temperature profiles in Fig. 8 are notably curved and follow a straight line only in the middle part of the wire.Although this effect is also present in Fig. 4, it is clearly much more pronounced here.However, many molecular dynamics studies have also found curved profiles near thermostatted regions (e.g.Refs.[6][7][8][9][10][11]).We believe that the extent of the curvature is determined by the mean free path of the heat carriers.As we pointed out in the introduction, some accommodation is expected near a thermostatted region and the longer the mean free path of the heat carriers, the longer such accommodations should reach into the wire.It is known that phonons with mean free paths above one µm contribute significantly to the heat transport in Si [24,25].This is consistent with the extent of the curvature shown by both models in Fig. 8. Due to the curvature of the temperature profiles, it is not possible to fit the complete profiles obtained from phonon Monte Carlo simulations to a straight line in order to derive the slope in the center.The green line in Fig. 8 was instead obtained by fitting a portion with length L/2 = 1 µm at the center of the wire.The same procedure was used to determine the profile slopes for other wire lengths.
In Fig. 9 we present the inverse slopes obtained from phonon Monte Carlo simulations with the two scattering models.In the case of Si-1 the inverse slopes follow a line parallel to the prediction by Fourier's law for wires longer than about 20µm.For shorter wires, the slopes begin to diverge from this line.Similarly, the inverse slopes obtained with model Si-2 follow a straight line for wires longer than about 2µm and diverge from this line for shorter wires.Qualitatively, this is the same behavior shown by the molecular dynamics simulations (cf.Fig. 6).As discussed above, we believe that the divergence from the straight line for shorter wires is a result of the fact that the profiles are not simply straight lines with jumps at the interfaces.As shown by the curvatures in Fig. 9 the temperature profiles reach a constant gradient only at a certain distance into the wire.For shorter wires the accommodations from both ends overlap resulting in a higher slope.

Thermal resistances
The phonon Monte Carlo method also allows the calculation of the heat flux in the wires which we use again as an independent test quantity.Fig. 10 shows the thermal resistances obtained from the heat fluxes in our simulations of Si nanowires.The figure shows that for both scattering models the resistances follow a straight line with a nonzero offset.Although this offset appears small in the figures, it adds a significant amount of resistance to shorter wires.This behavior is similar to the results obtained from molecular dynamics simulation (cf.Fig. 7).It is interesting to note that both simulation methods do not show noticeable deviations from a straight line for shorter wires.This means that the resistance added by the interface with the thermostatted regions remains unaffected by the mechanism leading to changes in the thermal gradients in shorter wires.
The values of the attenuation coefficient µ obtained from the offset of the fit lines in Fig. 10 are 1.05 µm −1 for model Si-1 and 3.69 µm −1 for Si-2.These values are substantially smaller than   The blue line is a linear fit of the data.the values 11.4 µm −1 and 161 µm −1 obtained from the fit lines in Fig. 9 for the two models.This leads to the conclusion that the radiator model provides a qualitatively correct picture but is not able to predict the value of the constant resistance offset based on a single attenuation constant alone.
A quantitative prediction of the offset probably needs to take into account the details of the mechanisms leading to the curved temperature profiles at the ends of the wires.

Summary and conclusions
In this work we have carried out non-equilibrium molecular dynamics simulations and phonon Monte Carlo simulations of thermal transport through nanowires.The results from these simulations are used to study the impact of ballistic transport effects on the temperature profiles and the thermal resistance of nanowires and similar quasi one-dimensional nanostructures.We also present a simple phonon radiator model.This model provides insight into the behavior observed in the simulations and an understanding of the underlying physical mechanisms.
In accordance with previous studies [5][6][7][8][9][10][11], our computer simulations show temperature gradients in the centre of the system that are reduced compared to the prediction of Fourier's law.In order to achieve the full temperature difference over the system length the temperature profiles are curved with steeper gradients at the end of the system close to the thermalized leads.The degree of curvature of the profiles varies between the simulations.We attribute the variation of the curvature in the simulations to differences in the phonon mean free path.It seems reasonable that in a system with a longer mean free path the effects of the thermal leads reach deeper into the system.
The radiator model predicts exactly linear temperature profiles with discrete jumps at the end of the system.While this is not exactly the behavior shown by the simulations, we believe that the steep gradients and curved profiles in the simulations are a continuous manifestation of the discrete jumps.Both simulation methods account for a range of phonon modes with different mean free paths whereas the radiator model only accounts for a single carrier particle with a uniform attenuation constant.Although the model does not exactly predict the behavior in the simulations, it provides a physical explanation for the reduction of the temperature gradients.Without this reduction the thermal flow in the wire does not match the flows generated by the thermal leads.
The inverse temperature gradients (profile slopes) in the centre of the wires obtained from our simulations follow a straight line for longer wires.This behavior is predicted by the phonon radiator model.For shorter wires the inverse gradients diverge from the straight line.A likely explanation for this are the curved tails of the temperature profile which reach a finite distance into the wires.For shorter wires the tails from both ends will begin to overlap resulting in an increased temperature gradient.
In addition to the temperature gradients, we determined the thermal resistance of the systems from our simulations.As predicted by the radiator model, the results follow a straight line with a non-zero offset for the full range of system lengths.The offset represents a constant thermal resistance caused by the interface to the thermal leads.While the variation of the thermal resistance is qualitatively predicted by the radiator model, it appears that the model is not able to predict the size of the resistance offset.Similar to the temperature gradients, a quantitative prediction would probably require a more sophisticated model.
The addition of a constant resistance to the wires by the thermal leads is similar to the wellknown Kapitza effect which causes a constant thermal resistance at the interface between materials.One might ask whether the artificial dynamics of the thermostats in the molecular dynamics simulations modify the vibrational properties in the thermostatted regions sufficiently to cause a Kapitza resistance.However, if this was true, why do the phonon Monte Carlo simulations (as well as the radiator model) show the same effect?In the phonon Monte Carlo simulations there is no difference between the materials in the thermal leads and the wire.Phonons enter the wire from the thermal leads with the expected frequency distribution and without any scattering so that there is no reason for a Kapitza resistance.Furthermore, our tests with different thermostat settings show that the temperature profile in the wire is unaffected by the details of the thermostat dynamics.This also speaks against the Kapitza mechanism as a cause for the additional resistance.
In summary, our simulations show the impact of ballistic transport on the heat conduction in nanowires.With some limitations, the radiator model is able to account for these deviations.For wires much longer than the typical mean free path of phonons, the model predicts merely a constant offset of the resistance which becomes negligible for wires of macroscopic length.For wires that are closer to the Casimir limit [26], i.e. wires that are much shorter than the typical mean free path of phonons, a more refined model that captures details of the scattering mechanism is required.
In our opinion, the ballistic transport effects studied in this work might be relevant for two reasons.The first reason is the additional term that appears in the thermal resistance.The second reason is the non-constant temperature gradient which might be important in applications that are sensitive to the local value of the temperature gradient.Both the steep gradient close to the thermal leads as well as the reduced gradient elsewhere might have to be considered in such applications.

Fig. 1
Fig. 1 Nanowire between two thermal leads.The upper part of the figure shows the temperature profile predicted by Fourier's law.The insets indicate the difference in the local environments in the middle of the wire and close to the cold lead.

Fig. 2
Fig. 2 Temperature profiles predicted by the radiator model for a wire of length L and various attenuation constants µ.The blue and red shaded areas indicate the position of the cold and hot thermal leads.

Fig. 3
Fig. 3 Simulation cell used in the molecular dynamics simulations with periodic boundary conditions.The upper part of the figure shows the temperature profile predicted by Fourier's law.

Fig. 4
Fig.4Average thermal profile in a wire of length L = 100 a 0 .The orange, blue and red dots indicate the temperatures in the wire, cold lead and hot lead sections of the system, respectively.The green line is the result of a linear fit of the data in the wire section.The thin black line shows the behaviour expected from Fourier's law.The blue and red shaded areas indicate the position of the cold and hot thermal leads.

Fig. 5
Fig. 5 Left half of the temperature profile in a wire of length L = 100 a 0 for different thermostat settings.

Fig. 6
Fig. 6 Slope (left) and inverse slope (right) of the temperature profile in the wire section as a function of the wire length L. Black dots indicate the results from molecular dynamics simulations.Orange lines show the results expected from Fourier's law.The blue lines are the result of a fit of the inverse slope data with L ≥ 200 a 0 .

Fig. 7
Fig. 7 Thermal resistance of the wire configurations as a function of the wire length.The blue line is a linear fit of the data.

Fig. 8
Fig.8Average thermal profile in a Si wire of length L = 2 µm from phonon Monte Carlo simulations using the scattering models Si-1 (left) and Si-2 (right).Temperatures in the wire, the hot lead and the cold lead are shown in orange, red and blue, respectively.The green lines are result of linear fits of the data in the center half of the wires.The thin black lines shows the behaviour expected from Fourier's law.The blue and red shaded areas indicate the position of the cold and hot thermal leads.

Inverse Slope m - 1 [μm K - 1 ]Fig. 9
Fig.9Inverse slope of the temperature profile in Si wires as a function of the wire length L obtained from phonon Monte Carlo simulations using the scattering models Si-1 (left) and Si-2 (right).Black dots indicate the simulation results.The orange lines shows the behaviour expected from Fourier's law.The blue line is the result of a fit of the inverse slope data for wires with L ≥ 20 µm (left) and L ≥ 2 µm (right).

Fig. 10
Fig.10Thermal resistance of Si nanowires as a function of the wire length obtained from phonon Monte Carlo simulations using the scattering models Si-1 (left) and Si-2 (right).