Molecular Dynamics Simulation of Heat Conduction through a Molecular Chain

Molecular Dynamics Simulation of Heat Conduction through a Molecular Chain

(Parte 2 de 4)

This might explain why the time constant τ2 is smaller for the species with n ) 0. Experimentally,4 a slow IVR component was not detected. However, the absolute value of τ2 is in the same range as typical VET times, so that, even if present in the real molecular systems, it might be difficult to detect. Therefore, in the following discussion, only the faster component is considered and used for comparison with experiment.

In addition to the alkyl chains, a methoxymethyl chain and 1,2-dimethoxyethylchain,as well as a thiomethoxymethylchain, were analyzed.The substitutionof a CH2 unit with oxygenslows the IVR and increases the corresponding time constant by roughly 40%. This is consistent with the experimental data for the propyl and dimethylether compound.18,23 Doubling the mass of the second unit in the chain by exchanging oxygen through sulfur further increases the IVR time constant by 40%. In one-

Figure 2. Ensemble-average temperature equilibration of azu-(CH2)3-ant dissolved in Xe. The temperature Ti(t) was averaged over all atoms belonging to the subunit i. The dashed, dotted, and solid lines represent i ) azu, alk, and ant, respectively.

14042 J. Phys. Chem. A, Vol. 113, No. 51, 2009 Schroder et al.

dimensional chain models, the mass effect was observed to be much more pronounced.25 However, in our results, the interac- tion potentials of adjacent CH2 units seem to attenuate the mass effect, and the deceleration of IVR seems to be caused by the smaller number of local modes transporting the energy through the chain and the reduced amplitudes of the local modes at the chromophore-chain interfaces.39

In Figure 4, the determinedIVR time constants,τ1, are plotted as a function of the alkyl chain length n and compared with experimental data. For long chains, the IVR rates derived from the simulations are roughly 50% slower than the experimental values. Nevertheless, the overall trend in the chain length dependence of τ1, at least for the azulene component, is well reproduced by showing a saturation of τ1 at alkyl bridge lengths with n > 4. For the anthracene moiety, this effect is less pronounced. We merely observe that the rise of τ1 with n slows as the chain length continues to increase.

Transient and Steady-State Temperature Profiles. The observed plateaus in the experimental and simulated τ1 values of azulenein Figure 4 suggestthat they are dominatedby similar energy-transfer mechanisms. The latter is analyzed here by studying the temperature profiles in the molecular chains during relaxation. This is exemplified in Figure 5 by means of the compound azu-(CH2)19-ant. The first and last circles in each graph refer to the azulene and anthracene temperatures, respec- tively. The circles between them belong to successive CH2 units in the alkyl chain. The profiles were obtained 0.5 ps (black circles), 5 ps (gray circles), and 20 ps (white circles) after azulene excitation. The simulation data show that, during the whole relaxation process, a strong temperature gradient at the azu-alk interface persists, whereas in the chain and at the alk-ant interface, this gradient vanishes. Only for shorter chains (not shown here) are strong temperature gradients observed also at the alk-ant interface, indicating that the coupling of the chromophores to the alkane chain is limiting the energy transport. In the chain, the flux is obviously much faster. The strong temperature gradient at the azu-alk interface in Figure 5 extends to the fourth methylene group. This explains why the azulene relaxation time in Figure 4 reaches a plateau at n > 4.

Analogous nonlinear temperature profiles have often been found in computer simulations30,31,45,46 of thermal conducting solids, for example, across grain boundaries in diamonds,29 and are closely related to the so-called Kapitza resistance.27,28 To investigate this effect in more detail and also determine a

Figure 3. Ensemble-average temperature equilibration for isolated systems azu-(CH2)n-ant with varying n. The temperatures Ti(t) of the units i ) azu, alk, and ant are illustrated by dashed, dotted, and solid lines, respectively.

TABLE 1: Characteristic Intramolecular Energy-Transfer Times in Bridged Azulene-Anthracene Compoundsa,b azu ant a Time constants for the azulene (azu) and anthracene (ant) sides were obtained from exponential fits. b Values in parentheses are experimental data.4,18

Heat Conduction through a Molecular Chain by MD J. Phys. Chem. A, Vol. 113, No. 51, 2009 14043

quantitative thermal transport parameter of the chain, steadystate simulations with thermostatted chromophores turn out to be more advantageousthan calculationswith transientexcitation as presented so far. For this purpose, the temperatures of azulene and anthracene

time of 0.1 ps. These values correspond to the chromophore temperatures of the simulations with transient excitation after harmonic dephasing before IVR starts. The steady-state temperature profiles are shown in Figure 6 for various values of n. In contrast to Figure 5, even for the longer bridges, temperature gradients in the chain and, in particular, at the alk-ant interface appear. We show below that the temperature gradient at the azu-alk Kapitza region is slightlylarger than that at the alk-ant interface. Obviously, when a steady-state flux is not enforced, as is the case for the simulations with transient excitation, only the interface with the strongest gradient limits the energy flux and appears in the temperature profile (see Figure 5).

Because the steady-state heat flux J is independent of time and position, the spatially dependent heat resistance is easily derived by means of eq 5. The temperature profiles of Figure

6 suggest that the CH2 chain be considered as being composed of three parts: two regions close to the chromophores, where the Kapitza effect is pronounced, and a central part, where the heat resistance is constant and equal to rc, the resistance of the homogeneous chain far from its ends. The central part is absent for short chains, for which the Kapitza regions are superimposed on each other. Theoretical estimations 27,47 suggest that the Kapitza resistance decreases exponentially with the distance from the interphase. Thus, we model the heat resistance by where ra, λa and rb, λb are constants corresponding to azu and ant, respectively, and L ) (n + 1)l is the total chain length.

We require the heat resistance to be a continuous function of the coordinate x. Then, the temperature is a continuously differentiable function of x, and it is a simple task to find the boundaries, xa and xb, of the region with constant resistance rc. They are defined as the points of the Kapitza regions where the spatial derivatives of temperature are equal to the constant temperature gradient in the central region.

The temperatureprofileis calculatedby integratingthe Fourier law at constantflux J. The results can be written for the three regions as

Figure 4. Energy loss and gain times for the azulene (full circles) and anthracene (gray diamonds) sides vs alkyl chain length (cf. Table 1). The squares represent the experimental values.18

Figure 5. Transient temperature profiles of the alkyl chain of azu-(CH2)19-ant including azulene (left) and anthracene (right) after 0.5 ps (full circles), 5 ps (gray circles), and 20 ps (white circles).

Monoexponential fits (dashed lines) provide the characteristic width of the Kapitza region at the azu-alk interface of 1.2 C-C bond length.

14044 J. Phys. Chem. A, Vol. 113, No. 51, 2009 Schroder et al.

where Ta and Tb are the temperatures of the boundaries of the constant-resistance region and

It is worth to notingthat the S-shapedsteady-statetemperature profiles cannot be explained by a temperature dependence of the heat conductivity. In this case, the derivative of eq 5 yields

The second derivative of the temperature has opposite signs near the chromophores. Thus, to satisfy eq 12, the temperature derivative of the thermal conductivity has to change its sign on the way from one chromophore to the other, which seems to be unrealistic.

The heat conductionequations(eqs 4, 5, and 7-9) are written for a continuous system, whereas we deal with a discrete chain model. Its grid points represent heat sources and reservoirs, and the bonds can be considered as heat resistances. The bond heat resistance is equal to the integral of r(x) over the bond length.

In the region of constant resistance, the bond resistance is Rc ) lrc. The Kapitza bond resistances near the azulene take the form

Because the Kapitza region extends up to five CH2 units, a meaningful value of the bond resistance of the central chain can be evaluated only for alkyl chains with n > 12. The temperature profiles of Figure 6 for these long bridges show well-separated central segments with constant temperature gradients of -1.3 ( 0.1 K/l. The heat flux J was computed by means of eq 6 between every methylene unit during the steadystate simulations and averaged over all interfaces. Consistency of the data analysis was verified by comparing the different interfaces with respect to constant J. In Figure 7, the heat flux is plotted versus the total number of CH2 units. With increasing chain length, the heat flux decreases up to a chain length of seven CH2 units. Beyond this length, the heat flux is almost constant and reaches an average value to 21.6 nW. The qualitative trend of J(n) is in accord with the dependence of the relaxation time on chain length found after transient excitation of azulene (see Figure 4). With J and the temperature gradients at hand, the central chain bond resistance and the total chain resistance can readily be calculated as Rc ) 1.3 K/J )

0.06 K/nW and R ) (Tazu - Tant)/J ) 20 K/nW, respectively. For long chains, both quantities are independent of n.

Exponential fittings to the long-chain Kapitza regions yield the chain lengths: For n ) 19, 16, and 13, respectively, we find asymmetry between the two sides of the chain is evident. It can be explained by the different natures and mode spectra of the two chromophores. A second reason might be the direction of the heat flux. To separate these two effects, steady-state simulations with identical chromophores on the two sides are in progress. The fitting constants are used to calculate the heat resistanceof those bonds connectingthe two chromophoreswith the chain. From eq 13, for n ) 19, 16, and 13, respectively, we obtain values of 6.9, 6.5, and 6.2 K/nW at the azulene side and 5.2, 5.7, and 6.2 K/nW at the anthracene side. These numbers are 2 orders of magnitude larger than Rc in the center of the chain.

The thermal diffusivity of the chain in the region of constant can be compared with the thermal diffusivity of highly oriented polyethylene parallel to the draw direction, for which Rc ) 9 × 10-6 m2/s was measured.48 The reasonfor the differencemight be that, in classicalsimulations,all vibrationalmodes are excited and contribute equally to the heat capacity whereas, in the polyethylene chain, high-frequency modes are not excited.

Figure 6. Steady-state temperature profiles of azu-(CH2)19-ant with n ) 3, 5, 7, 9, 13, 16, and 19. The temperatures of the two chromophores were stabilized at 750 and 320 K, respectively, using Nose-Hoover thermostats. The first circle (x ) 0) corresponds to the azulene temperature; the last circle of each curve corresponds to the anthracene temperature.

Heat Conduction through a Molecular Chain by MD J. Phys. Chem. A, Vol. 113, No. 51, 2009 14045

Hence, the average heat capacity cL of a methylene unit is significantlysmallerthanassumed,suchthatRc is underestimated. Molecular dynamics simulations of heat conduction in solids31,46,49-51 have shown strong Kapitza effects at the interface betweenthermostattedand nonthermostatedregions.In principle, such artificial Kapitza effects could also be present in our simulations when the chromophores were thermostatted to establish a steady-state heat flux. The importance of this effect can be estimated by comparing the heat flux values of these simulations with those of our time-variant simulations, where, in the course of relaxation, no thermostats were used. The latter can be calculated according to J ) |dEazu/dt| = 51(kB/τ)∆Tazu(t), where the 51 degrees of freedom of azulene are considered and

∆Tazu = 250 K is the temperature difference between azulene and the alkyl chain in the steady-state simulations. With τ = 8 ps for the long chains (cf. Figure 4), the heat flux is calculated to be approximately 2 nW, which is very close to the value of 21.6 nW derived from the steady-state simulations. Therefore, we can conclude that, in our case, artificial Kapitza effects are of minor importance.

For heat transfer in strictly one-dimensional chains, another type of anomalywas reported.25 In computerstudies,the thermal conductivity of long but finite chains was found to diverge as n , where the coefficient varies between 0.35 and 0.4. In our case, however, the chain is embedded in three-dimensional space. This might be the reason that the heat diffusivity and heat conduction coefficient for the long chains are constant and this effect is not observed.

Summary and Conclusions

In this study, vibrational energy transfer between the two chromophores azulene (azu) and anthracene (ant) connected by a molecularchainof variablelengthwas investigatedby classical MD simulations. The chain involves 0-19 methylene units; likewise,ether and thioethercompoundswere examinedas well. In a first series of simulations, the relaxation toward intramo- lecular equilibrium after selective excitation of the azulene side was investigated. Analogous to a recent laser experiment, an excess energy of 18000 cm-1 was statisticallydistributedamong the vibrationalmodes of the azuleneunit. The subsequentenergy flux through the bridge to the anthracene side was followed by monitoring of the temperatures of the two chromophores and each methylene unit. Comparison of the relaxation of azu-(CH2)3-ant in liquid xenon and under isolated conditions reveals that the influence of the solvent in the intramolecular energy flow is weak. Therefore, further analysis concentrated on simulations of these molecules in the gas phase.

For most of the molecules,the temperaturedrop at the azulene side could be modeled exponentially. The corresponding time constants are in the range of a few picoseconds and show a similar dependence on the number of methylene units n of the alkyl chain as in the corresponding laser experiments: For short bridges, the relaxation time increases proportionally to n, whereas for n > 4, it levels off and becomes independent of n. A closer analysis of the temperature profiles in these molecules indicates a pronounced Kapitza effect at the azulene-chain interface characterized by strong thermal resistances of the first few methylene units. As a result, the interface and not the remaining chain acts as a bottleneck for the energy drain of the excitedazulenesuch that the correspondingtime constantsabove n ) 4 are independent of chain length.

To quantify the C-C bond resistancesRc in the alkyl bridges, additional calculations with thermostatted chromophores were performed, providing steady-state temperature profiles with constant energy flux. It was found that, on both sides, Rc decreasesexponentiallywith the distancefrom the chromophore and levels off at the fifth methylene unit. For long chains (i.e., n > 12), a domain with constant Rc is observed in the center of the bridge. Here, the bond resistance is 2 orders of magnitude smaller than at the edges of the Kapitza regions. Contrary to true one-dimensional chains, the heat conductivity coefficient in the bridge center does not depend on chain length, probably

Figure 7. Average steady-state heat flux during the simulation with thermostatted chromophores as a function of alkyl chain length. 14046 J. Phys. Chem. A, Vol. 113, No. 51, 2009 Schroder et al.

because our model systems are embedded in three-dimensional space. The resulting heat diffusivity compares reasonably well with that of highly oriented polyethylene in the draw direction. It is shown that the observed temperature profiles cannot be explained by a temperature dependence of the heat conductivity and that artificialKapitzaeffectsat the thermostat-chromophore interface are of minor importance.

Acknowledgment. C.S. thanks Prof. O. Steinhauser and Dr.

A. Svrcek-Seilerfor very helpful discussions.The authors thank one of the reviewers for constructive criticism and suggestions. Financial support of this project by the Deutsche Forschungsgemeinschaft (SFB “Molekulare Mechanismen unimolekularer Prozesse”) is gratefully acknowledged.

Appendix A Parameterization.

Parameterization of the force field in terms of bond, angle, and dihedral potentials was first done separately for the two chromophores (azu, ant) and the saturated chain (alk), taking advantageof the smallernumberof normalmodes to be adjusted for the subunits, where azulene, anthracene, and butane served as the respective model bases. In a second step, the interfacial bonding potentials at the azu-alk and alk-ant boundaries were determined so as to reproduce the theoretical normal-mode spectrum of the bridged azu-alk-ant molecule. The normalmode vibrational frequencies of subunits and bridged molecules were computed by the B3LYP method within DFT using the Gaussian 98 program package52 and scaled by the empirical factor 0.9613.53 It turned out that neither the intrachromophore and intrachain potentials nor the interfacial force field needed to be changed in order to reasonably reproduce the theoretical mode spectra.This has the additionaladvantagethat a systematic comparison of energy redistribution in bridged bichromophoric molecules of varying alkyl chain length is possible without having to consider possible effects of changes in the force-field parameters on the time scales and mechanism of IVR.

Atomic partial charges as obtained from the Gaussian 98 program output were not used in the molecular dynamics simulations, because (long-range) polar interactions were not regarded as essential for IVR in the isolated molecules studied here and for intermolecular energy transfer in xenon solvent. The inclusion of intramolecular van der Waals interactions did not change the simulation results significantly for a subset of samples with varying chain length. Therefore, in the systematic variation of alkyl chain length, these types of potentials were omitted in order to additionally enhance the simulation speed.

However,we also simulatedthe azu-(CH2)3-ant system in 225 Xe atoms in order to analyze the interplay of IVR and VET in our systems. In this case, we used the Lennard-Jonesparameters σ ) 4.055 Å and ε ) 0.457 kcal/mol for the interaction of the carbon atoms with Xe. Furthermore, these nonbonded interactions complicate the parametrization with respect to the normal modes because strong intermolecular potentials shift the normal modes to lower frequencies (even imaginary frequencies were detected). In addition to the parametrization, the interpretation of the IVR is much more complicated. Because of the nonbondedinteractions,parts of the excessenergyof the hot azulene are transferredto the anthracene without taking the path through the alkyl bridge. Another reason concerns the analysis of the heat current. In future equilibrium simulations, we plan to evaluate the heat flux through the chain directly. A detailed Green-Kubo analysis seems far to complex if each chain atom interacts with all other atoms in the system.

TABLE 2: Force Field of Azulenea Bonds



(Parte 2 de 4)