Thermal transport in fullerene-based molecular junctions: molecular dynamics simulations

We investigate phonon thermal transport of fullerene-based single-molecule junctions by employing classical molecular dynamics (MD) simulations. We compute the thermal conductances of C60 fullerene monomers, dimers, and trimers utilizing three distinct MD methods. We observe the equilibration dynamics in one approach, and employ two other nonequilibrium steady state simulation methods. We discuss technical aspects of each simulation technique, and show that their predictions for the thermal conductance agree. Our simulations reveal that while the thermal conductance of fullerene monomer and dimer junctions remains similar, that of trimer junctions experiences a significant reduction. This study could assist in the design of high-performing thermoelectric junctions, where low thermal conductance is desired.


I. INTRODUCTION
Recent groundbreaking experiments reported the thermal conductance of selected single molecules in the configuration of a metal-molecule-metal junction [1][2][3] .These efforts are of an interest from the fundamental side for understanding thermal transport mechanisms at the nanoscale level of a single molecule [4][5][6][7] , and for applications in electronics, energy conversion, and novel materials development [8][9][10][11][12] .In particular, for flexible organic-based electronic devices, it is important to identify molecules that exhibit large thermal conductance, beneficial for transferring away dissipated heat and maintaining the integrity of the junction.Conversely, thermoelectric devices show improved efficiency when reducing their thermal transport, including the electronic, phononic and photonic components.
In fact, due to the challenge in measuring thermal transport of single-molecule junctions, simulations play a leading role in guiding the search for effective thermal transport junctions 29 .Two distinct computational tools dominate this field: (i) Employing a quantum method that integrates density functional theory (DFT) with the nonequilibrium Green's function (NEGF) technique 22,30 .This framework does not take explicit many-body interactions into account.It enables separate computations of the contributions of electrons or phonons to the thermal conductance.The DFT-NEGF method builds on the Landauer equation, and it assumes noninteracting electrons and harmonic interactions in the junction, ignoring dynamical effects: Simulations are executed on a frozen configuration, yet frequencies and metal-molecule hybridization energies are obtained from DFT calculations.(ii) Using classical molecular dynamics (MD) simulations that build on parametrized force fields.Conventional classical MD simulations cannot provide any information on electronic thermal transport.However, for the vibrational-phononic component, MD simulations offer great flexibility as such methods are not restricted to the harmonic force field.Furthermore, MD simulations naturally allow the molecule and the metal atoms to move throughout the simulation, thus taking into account dynamical effects.However, in the standard implementation of MD simulations, such as that used in the LAMMPS simulator 31 , potentials are phenomenological with fixed parametrization.In our study here we employ such an MD approach with details given in Sec.II.
Fullerenes, carbon-based cage-like molecules, support a wide range of applications in nanoscience in e.g., drug delivery, catalysis, organic photovoltaics, and energy storage devices 32 .In particular, electrical conductance and thermopower measurements of C60 monomer and dimers were carried out in Ref. 33.It was further suggested in that study that stacks of C60 molecules could be candidates for thermoelectric applications by achieving high efficiency.However, a DFT-NEGF computational study of thermal transport in monomer and dimer fullerenes did not support this proposal; it was shown in Ref. 34 that phonons dominate and limit the thermoelectric performance of dimer fullerene junctions.
Motivated by thermopower measurements of fullerene monomers and dimers metal-molecule-metal junctions, and the study of their electronic, thermal and thermoelectric prop-Figure 1. Visualization of a single-molecule gold-fullerene-gold nanojunction (OVITO 35 ) and the methods used to simulate thermal transport as described in Sec.II.(Top) In the AEMD method, the two metals are initially separately thermalized to T h and T c , while the molecule is prepared at the average temperature.The thermal conductance is evaluated from the equilibration dynamics, see Sec.II A. (Centre) In the LNEMD method, Langevin thermostats are attached to a certain region at the back of the metals away from the junction, thermalizing them to temperatures T h (left) and T c (right).The heat current is evaluated by calculating the net heat exchange between the thermostats and the metal at each contact, see Sec.II B. (Bottom) In the RNEMD method, the heat current into the system (and equivalently out) is imposed as a parameter, resulting in a temperature drop on the junction.The method is described in Sec.II C. In all three approaches, a group of 270 gold atoms at the left and right boundaries (colored in yellow) are fixed, to keep the structure from collapsing.erties using the DFT-NEGF approach, we simulate here the phononic thermal conductance of fullerene junctions using the complementary classical MD approach.Contrary to the Landauer equation, classical MD allows us to include anharmonic interactions and dynamical effects.Our first main objective in this study is the computation of the thermal conductance of fullerenes in a junction configuration, including monomers, dimers and trimers junctions.Concurrently, our second goal is to continue establishing molecular dynamics simulation methods for the study of single-molecule thermal transport.We achieve this goal by adopting and comparing three molecular dynamics approaches for thermal transport calculations.In this respect, while in our current calculations we utilize parameterized force field potentials and rely on fully classical simulations, more advanced studies could benefit from the adoption of first-principle machine learning potentials and the integration of quantum corrections.The confidence established in simulation methods through this study will serve future investigations that aim to employ more accurate potentials and incorporate quantum corrections.
The plan of this paper is as follows.In Sec.II, we describe three molecular dynamics methods for the simulation of thermal transport in single molecule junctions: The Approach to equilibrium MD, which follows equilibration dynamics, the Langevin nonequilibrium MD, in which we thermostat the end sections of the two metals to different temperatures, and compute the net heat current in the system, and the reversed nonequilibrium MD technique, where we do the opposite: We impose the net heat current and compute the temperature profile.Simulation results are presented in Sec.III, where we discuss the heat current trend with temperature difference and the number of stacked fullerenes.We summarize our work in Sec.IV offering future perspectives.

II. SETUP AND SIMULATION TECHNIQUES
We construct metal-molecule-metal junctions and perform thermal transport simulations using three methods: Approach to Equilibrium Molecular Dynamics (AEMD), Langevin Non Equilibrium Molecular Dynamics (LNEMD), and the Reversed Non Equilibrium Molecular Dynamics (RNEMD), as illustrated in Fig. 1.We calculate the thermal conductance of the junction G in the three methods.In addition, the latter two nonequilibrium methods enable us to study the heat current in the steady state limit as a function of factors such as the imposed temperature difference.
The junction consists of two gold leads, each containing 2160 atoms.At the ends of each lead, the Au atoms are fixed in order to hold the junction and prevent the two pieces of gold from collapsing together.These atoms cannot move in the simulation.Periodic boundaries are implemented in the x and y dimensions (planes of gold), but the simulation box is extended in the z dimension; the top of the box is placed farther away from the atoms.
Atom types comprise of Au in the gold leads and C in the fullerene molecule placed between the two leads, see Fig. 1.The force fields used are the same for all three methods.Interactions within the fullerene molecule are modelled using an adaptive intermolecular reactive bond order (AIREBO) potential that is optimized for carbon-carbon interactions 36 .Au-Au interactions are given by the embedded atom model.A Lennard-Jones potential is set for long-range, non-bonded interactions between all atoms.The parameters for the potentials are summarized in Table I.

A. Approach to equilibrium MD
The approach to equilibrium molecular dynamics was described in several studies, for example Refs.27,39,40.In this method, the thermal conductance is calculated from the equilibration dynamics by monitoring the temperatures of the two bulk solids.The method relies on the phenomenological Newton's law of cooling.We start with a nonequilibrium condition where the two separate metals are prepared at distinct temperatures.Once the metals are attached via the fullerene at the junction, energy is exchanged between the metals and the system eventually reaches a global equilibrium state.By monitoring the temperature of the two bulk metals as they equilibrate, one can determine the thermal conductance of the junction as follows: Newton's cooling law for each metal is given by Here, the temperatures change in time towards equilibrium.C is the heat capacity of the the metal.Using the Dulong-Petit law, it is approximated as C = 3k B N Au with N the number of moving gold atoms in that metal; k B is the Boltzmann's constant.Defining the temperature difference, ∆T (t) = T h (t) − T c (t), the two equations ( 1) can be combined, The solution to this first-order differential equation is where τ ≡ C/2G is the equilibration timescale.We also define the averaged temperature as T (t) = [T h (t) + T c (t)]/2.This temperature is evaluated as the average over the full bulk.

Procedure
In the AEMD method, NPT equilibration is first carried out to relax the system and the simulation box, barostatting for zero pressure and thermostatting to a target temperature, T .Subsequent NVT equilibration is carried out separately on the two Au leads to achieve the starting temperatures, T h (0) and T c (0).The fullerene molecules at the center are equilibrated to the averaged value T .These equilibrations are performed over 1.5 ns with 1 fs timesteps.As for the production run: An NVE simulation is conducted for 5 ns, so the temperatures of the Au relax to the equilibrium value.Given the symmetry of the setup, the temperature of each metal, and the fullerenes, eventually reaches T = [T h (0) + T c (0)]/2.Temperature data is logged every 1000 steps of 1 fs timesteps, equivalent to every 1 ps.

Data Analysis
We display an example of AEMD simulations in Fig. 2 with initial temperatures of 125 K and 75 K at the two metals.In panel (a), we follow the equilibration dynamics, which takes place within about 2 ns after contact.After that time, the temperatures of the two metals reach the averaged value (100 K), showing some temporal fluctuations.Temperatures are evaluated from the kinetic energy assuming equipartition.The temperature of each metal is calculated as the average of the bulk gold (1890 atoms).
To extract the relaxation timescale as described in Eq. ( 3), we plot in Fig. 2 (b) the evolving temperature difference as a function of time, illustrating an exponential decay within the first ≈2 ns.In fact, we only use the data marked in black (up to 1.5 ns) for the calculation of the thermal conductance, with results summarized in Fig. 4. We repeat this calculation at two different averaged temperatures, T = 100 K and T = 300 K, with ∆T = 50 K and ∆T = 100 K respectively, in both cases reaching G ≈ 50 pW/K.
The primary advantage of this method lies in its relatively short simulation times.However, the relaxation dynamics may deviate from Newton's cooling law, displaying multiple relaxation times (as illustrated below in Fig. 7).While offering insights into nanoscale dynamics, these nontrivial trends pose challenges in calculating the thermal conductance of the junction.

B. Langevin nonequilibrium MD
In the Langevin Non-Equilibrium MD approach, several layers of gold atoms at each metal are continuously thermostated to a given temperature, which is different at the two ends, thus generating a nonequilibrium steady state, see Fig. 1(b).Examples of recent studies adopting this method for thermal transport calculations in polymer junctions include Refs.24-26.We use Langevin thermostats characterized by their temperature, T h and T c , with the damping timescale γ −1 .Here, γ corresponds to the rate of thermalization with the local Langevin bath.The instantaneous heat current at each contact can be evaluated from the Langevin dynamics, by building the trajectory ∆E cold (t) and ∆E hot (t) as the cumulative heat exchanged  3).The linear fit is performed over data up to 1.5 ns (black), while beyond that (gray) the system reaches thermal equilibrium.
from time t = 0 to time t.Our sign convention is such that heat input from the system's atoms to the attached thermostats is defined as positive, so ∆E cold (t) is positive while ∆E hot (t) is negative.

Procedure
In the LNEMD method, a temperature bias is imposed on the two gold leads using Langevin baths.NPT equilibration is first carried out to relax the system and simulation box, barostatting for zero pressure and thermostatting to a target temperature, T .A subsequent NVT equilibration is carried out by equilibrating the whole system to the same target temperature.Each equilibration is performed over 1.5 ns with 1 fs timesteps.As for the production step, NVE simulations are performed to produce steady-state temperature profiles.At this stage, Langevin thermostats are applied to the two gold leads, one set at a high temperature T h and the other at a cold temperature T c .The damping rate γ −1 is set to 0.05 ps, which was demonstrated in Ref. 28 to minimally interfere with the junction's dynamics.These production runs are executed for a total of 10 ns, with results averaged over the last 9 ns.Temper-ature data is logged every 1000 steps of 1 fs timesteps, equivalent to every 1 ps.

Data Analysis
Fig. 3 exemplifies data collection and analysis in the LNEMD method.Raw data of heat exchange with the thermostats is displayed in Fig. 3(a).The amount of heat absorbed by the gold atoms attached to the hot baths (red) is equal to the amount of heat released to the cold thermostats (blue).In the nonequilibrium steady state, the junction responds with the development of a temperature bias internally on the junction.In Fig. 3(b) we display this temperature profile.The temperature of each carbon atom is evaluated from its kinetic energy, and as expected, this measure for temperature shows large fluctuations (see error bars).For the gold metal, we evaluate temperature for slabs with 270 atoms.The temperature drop ∆T Au-Au is further marked in Fig. 3(b), evaluated at the edge gold atoms.The heat current in the junction is calculated from the slope in panel (a) as J = (|∆E(t 2 )| − |∆E(t 1 )|)/(t 2 − t 1 ), with t 2 > t 1 and for a long time interval.The thermal conductance is given by the ratio G = J/∆T Au-Au .As can be seen in Fig. 3(b), the thermal resistance is dominated by the metalmolecule interface 41,42 .To enhance thermal conductance, this boundary effect could be alleviated by e.g., using other materials as bulk contacts.
The LNEMD method is convenient to use, specifically if one is interested in studying the heat current-temperature difference trends e.g., towards studying the thermal diode effect and in investigations of the heat current fluctuations 28 .One should be aware however of two issues: (i) Properties of the molecular junctions such as its thermal conductance should not depend on the damping rate γ, and this should be verified for each system.This simulation parameter corresponds to the Au-thermostat interaction, and it is not part of the atomic model.(ii) Standard simulators such as LAMMPS do not implement solvers for stochastic differential equations, leading to violations of e.g., heat exchange fluctuation relations 28 .

C. Reversed nonequilibrium MD
In the RNEMD method, we use a reversed procedure for determining the current-thermal bias relationship.In this approach, the magnitude of the heat current is imposed as a parameter with J h energy per unit time added on average to one side-and the same amount extracted from the opposite metal contact.When simulating for a long enough time, the system reaches a nonequilibrium steady state building a temperature bias in the junction in response to the input heat current.The method was used extensively in recent molecular dynamics studies, even when combined with quantum corrections.Some recent examples include Refs.27,43-46.

Procedure
In the RNEMD method, heat is inputted and outputted close to the ends of the two gold leads, see Fig. 1(c).NPT equilibration is first carried out to relax the system and the simulation box, barostatting for zero pressure and thermostatting to a target temperature, T .Subsequent NVT equilibration is carried out by equilibrating the whole system to the same target temperature.Each equilibration is performed over 1.5 ns with 1 fs timesteps.In the production step, NVE simulations are performed to produce steady-state temperature profiles.Heat is inputted to one lead's bath and extracted from the other (as kinetic energy per unit time to the atoms), at the same rates.These production runs are executed for a total of 22.5 ns, with results averaged over the last 20 ns.Temperature data is logged every 1000 steps of 1 fs timesteps, equivalent to every 1 ps.

Data Analysis
In the steady state limit, we observe that in response to the enacted heat current, a temperature bias develops on the junction, similar to the profile shown in Fig. 3(b) for the LNEMD method.We define the resulting temperature bias on the junction as ∆T Au-Au ; the temperatures are evaluated on the metal at the proximity of the molecule.This thermal bias allows the calculation of the thermal conductance G using the linear response expression, G = J ∆T .Altogether, in the "reversed" approach we identify the J-∆T (current-temperature bias) relationship by setting the current and extracting the thermodynamic force, opposite to what is done in the direct, LNEMD method.
The main advantage of the RNEMD is that simulations are conducted in the NVE ensemble, which is straightforward to simulate.However, when aiming to understand the currentbias characteristics, especially in the context of the diode effect, the method becomes less convenient.This additional complexity arises due to the need to pre-estimate the temperature drop corresponding to a given current 27 .

III. RESULTS: MONOMER, DIMER AND TRIMER FULLERENE JUNCTIONS
We adopt the three molecular dynamics approaches, the AEMD, LNEMD and RNEMD, and simulate the thermal conductance of fullerene monomer, dimer and trimer junctions.Overall, we find that the three methods agree in their predictions.Particularly, the RNEMD and the LNEMD are similar in computational effort and in their results.Since the LNEMD is more intuitive to operate (one sets up the temperature at the boundaries and generates a heat current in the system), we overall find it to be the most convenient tool among the three.We now describe our results.
In Fig. 4(a) we display the thermal conductance of a C60 monomer nanojunction using the three different techniques at two different averaged temperatures, 100 K and 300 K. We make the following observations: (i) The thermal conductance of fullerene monomer gold nanojunction does not depend on temperature in the examined range of T =100-300 K. (ii) The three methods agree, providing thermal conductance of G ≈ 50 pW/K (49.62 pW/K).(iii) Repeating simulations showed small errors of less than 3 pW/K.(iv) In Fig. 4(b), we furthermore study the behavior of the heat current as a function of the temperature difference, ∆T Au-Au , and observe a linear trend.
In Fig. 5 we display the thermal conductance of stacked fullerene junctions with N = 1, 2, 3 units.We present results using two methods, the AEMD and the LNEMD, and results agree within the uncertainty of simulations.Significantly, we find that the thermal conductance of fullerene monomer and dimer junctions are comparable, at approximately 50 pW/K and 40 pW/K, respectively.Fullerene trimers however exhibit significantly lower conductance, measured at 10 pW/K.
Next, in Fig. 6 we display the heat current as a function of the applied temperature difference using the LNEMD method.We also calculate the slope, which corresponds to the thermal conductance.We find that even far from equilibrium, for ∆T / T ≈ 1, stacked fullerene junctions display linear heat current-temperature bias characteristics.
As for the size of the gap d between the metals, as one can see in the illustration of Fig. 5, the fullerenes lie in close proximity to each other (efforts to space them out led to their relaxation to this distance); the values d 1 and d 2 for the monomer and dimers, respectively, are in accord with those tabulated in Ref. 34.
The choice of method in future studies depends on specific research objectives.The AEMD has the clear advantage with its shorter simulation time for extracting the thermal conductance.However, as we also discussed in Ref. 27, Eq. (3) relies on substantial assumptions, most importantly that a single exponential decay characterizes the equilibration dynamics.
To assess this assumption, we present in Fig. 7 equilibration traces of the AEMD method for dimer and trimer junctions, complementing Fig. 2. Interestingly, while the dimer fullerene seems to follow a single exponential decay, the equilibration process of a trimer nanojunction is clearly more involved, manifesting two timescales: A relatively fast equilibration (yet slower than in the monomer and dimer junctions), followed by a slower decay.We note that we extracted the thermal conductance of Fig. 5 from the initial fast decay, indicated in black in Figs.7(b) and (d).The AEMD method offers nontrivial information on the equilibration dynamics, yet for thermal conductance calculations we suggest the LNEMD and RNEMD steady-state methods.These methods do not make any assumptions on the underlying dynamics and are thus more general and robust.We now compare our results to other studies.In Ref. 34, the authors computed the phononic contribution to the thermal conductance quantum mechanically, using the DFT-NEGF method, which however is limited to harmonic interactions.Furthermore, in their calculations, the junction was held at a fixed configuration.This is to be contrasted by our fully classical-dynamical simulations that yet include anharmonic interactions and long-range Lennard-Jones interactions.Another noteworthy difference between Ref. 34 and our work is the geometry of the metal contact.While here we form a junction between two flat metal surfaces, the junction was previously assumed to be formed between two tips, which can be sharp or blunt.
In Table II, we present results from Ref. 34 focusing on junctions with two blunt tips.This study further demonstrated that for monomers, the conductance was halved if one of the contacts was made sharp.The Table also summarizes our findings from Fig. 6.Interestingly, both quantum 34 and our classical simulations of monomer fullerene yield similar values for the phononic thermal conductance, at ≈ 45 − 50 pW/K, with the classical MD bringing results that are about 10% higher than the quantum ones, a typical observation as classical simulations often overestimate thermal transport.
DFT-NEGF and MD results though disagree on dimer junctions.While DFT-NEGF calculations provided a low thermal conductance at 7.0 pW/K 34 , our simulations predict that dimers conduct with a similar efficacy to monomers with a thermal conductance that is about 20% smaller than the monomer, at ≈ 40 pW/K.The strong suppression of thermal transport shows up only in the trimer junction, where MD predicts thermal conductance at ≈ 10 pW/K (we do not have results from the literature to compare to this value).We analyzed the temperature profiles of dimer and trimer fullerenes, similar to the monomer data presented in Fig. 3.We observed (not shown) a very small temperature drop between the fullerenes, at the order of 0.5 K for the dimer (and smaller for the trimer).Thus, internal contact resistance between the fullerenes cannot explain the substantial reduction in thermal conductance for trimer junctions.
To address the question of why there is a large difference in thermal conductances, we plot in Fig. 8 the vibrational density of states (VDOS) for isolated fullerenes based on velocity data from LAMMPS 47,48 .Below 6 THz, which is roughly the Debye frequency of gold, we anticipated the monomer VDOS to exhibit more peaks than the trimer's.In contrast, the trimer's VDOS displays more structure at low frequencies, likely representing modes associated with inter fullerene motion.However, these modes do not appear to contribute to the metal-to-metal thermal transport.The suppression of thermal conductance for trimer junctions thus remains an open question.
Going back to Table II, we comment that other experimental and computational studies on the thermal transport properties of fullerenes have explored a variety of configurations, rather than the metal-molecule-metal junction considered in our study.The wide variety of studied systems include fullerenes in solution 39,49 , crystalline fullerenes 50 , fullerenes with molecules confined within them 51 , and fullerenes manufactured into materials 52,53 .Furthermore, these studies typically focused on bulk thermal conductivity rather than thermal conductance, reflecting the consideration of fullerenes as a carbon-based nanomaterial.In the above mentioned systems, reported thermal conductivity for C 60 ranged from 0.1-0.4W/m K 49,50,52-54 .

IV. DISCUSSION AND SUMMARY
In this study, we investigated phononic thermal transport in single-molecule gold-fullerene-gold nanojunctions based on classical molecular dynamics simulations.We computed the thermal conductance of monomer, dimer and trimer fullerene molecules sandwiched between metals.We showed that the thermal conductance decreased monotonically with the number of fullerene units.Notably, the thermal conductance of trimer junctions was five times smaller than that of monomer junctions, and four times smaller than dimer junctions.
Our study further affirmed the consistency of different classical MD tools in studying thermal transport properties of nanojunctions.Particularly, we showed that the AEMD method, which extracts the thermal conductance values from transient dynamics, recovered correct results, and even when the relaxation dynamics followed a multi-exponential decay form.The next stage of such simulations would be to employ first principle-based parametrization for the atomic potentials, and the inclusion of quantum effects through corrections to the thermostats, as done in Ref. 25.These extensions require understanding of the underlying MD method, thus consistency checks performed here serve as a foundation for such more sophisticated modelling.
There are a plethora of open questions and simulation geometries that could be explored in future work.For example, it was demonstrated in ref. 34 that the metal-to-metal separation strongly affects the heat current.In the present work, we could not examine this aspect as the bulk metal electrodes tended to approach each other upon equilibration.One could however enforce a certain metal-to-metal gap by fixing (making immobile) a large portion of atoms at the boundaries.Additionally, the structure of the tip plays a significant role in thermal transport as illustrated in Ref. 34.This aspect was not examined in our work and we studied only the case of fullerenes placed between flat surfaces (which could represent very blunt tips).Finally, we contained ourselves here with gold metals, but results should depend on the type of the contact, particularly when using graphene electrodes instead of gold 25 .Other interesting extensions concern thermal transport characteristics of endohedral fullerenes, using encloses atoms or molecules as knobs to control thermal transport 51 .It is our hope that this work will spur additional computational and theoretical studies of heat transport in single-molecule junctions, and motivate new experimental efforts in this field.

36 Figure 2 .
Figure 2. Example of AEMD simulations of thermal conductance of C60 monomer nanojunctions.(a) The two gold metals are prepared initially at two different temperatures T h (red) and T c (blue), and these temperatures are monitored over time once the metals are put into contact through the fullerene molecule.(b) Data analysis of an AEMD simulation.A plot of ln(∆T (t)) over time provides the decay timescale τ, which is used for the determination of the thermal conductance of the junction through Eq. (3).The linear fit is performed over data up to 1.5 ns (black), while beyond that (gray) the system reaches thermal equilibrium.

Figure 3 .
Figure 3. Data collected in the LNEMD method.(a) Cumulative energy exchange of each metal with the Langevin thermostats.∆E hot and ∆E cold are the cumulative heat absorbed and emitted by Au atoms connected to the hot and cold thermostats, respectively.∆E ave (t) is the average value of these two trajectories (maintaining the sign positive for heat entering the thermostats).(b) Temperature profile across the junction in the steady state limit.Sections correspond to 270 Au atoms and single C atoms in the fullerene.Simulation parameters are T h = 125 K, T c = 75 K, and γ −1 = 0.05 ps.

Figure 4 .
Figure 4. Thermal transport in monomer fullerene junctions.(a)Thermal conductance using different methods as indicated at T =100 K (blue) and T =300 K (orange).(b) Heat current as a function of temperature bias in the RNEMD and the LNEMD methods at T =300 K.We repeated simulations for ∆T = 100 K and found that uncertainties in the heat current were smaller than the size of the data point symbol.

Figure 5 .
Figure 5. Thermal conductance of monomer, dimer and trimer fullerene junctions as depicted in the illustrations.The length of the metal-to-metal gap is further indicated; deviations during simulations were less than 1%.We compare LNEMD simulations to AEMD results showing a good agreement.

Figure 6 .
Figure 6.Heat current in stacked fullerenes as a function of temperature difference with N = 1, 2, 3 fullerenes top to bottom.The average temperature is 100 K.The slope of each curve, corresponding to the thermal conductance, G N , is indicated in each case.

Figure 7 .
Figure 7. AEMD simulations of (a)-(b) dimer and (c)-(d) trimer fullerene nanojunctions.In both cases, we display the temperature of the metals as a function of time, and the corresponding analysis to extract the thermal conductance.

Figure 8 .
Figure 8. Vibrational density of states (VDOS) for monomer and trimer fullerene as a function of frequency.The average temperature is 300 K.The black dashed line roughly marks the Debye frequency of gold, while the dashed blue line stands at room temperature.

Table II .
Table of computed thermal conductance of fullerene junctions