Molecular heat transport across a time-periodic temperature gradient

The time-periodic modulation of a temperature gradient can alter the heat transport properties of a physical system. Oscillating thermal gradients give rise to behaviors such as modified thermal conductivity and controllable time-delayed energy storage that are not present in a system with static temperatures. Here, we examine how the heat transport properties of a molecular lattice model are affected by an oscillating temperature gradient. We use analytical analysis and molecular dynamics simulations to investigate the vibrational heat flow in a molecular lattice system consisting of a chain of particles connected to two heat baths at different temperatures, where the temperature difference between baths is oscillating in time. We derive expressions for heat currents in this system using a stochastic energetics framework and a nonequilibrium Green's function approach that is modified to treat the nonstationary average energy fluxes. We find that emergent energy storage, energy release, and thermal conductance mechanisms induced by the temperature oscillations can be controlled by varying the frequency, waveform, and amplitude of the oscillating gradient. The developed theoretical approach provides a general framework to describe how vibrational heat transmission through a molecular lattice is affected by temperature gradient oscillations.

The time-periodic modulation of a temperature gradient can alter the heat transport properties of a physical system.Oscillating thermal gradients give rise to behaviors such as modified thermal conductivity and controllable time-delayed energy storage that are not present in a system with static temperatures.Here, we examine how the heat transport properties of a molecular lattice model are affected by an oscillating temperature gradient.We use analytical analysis and molecular dynamics simulations to investigate the vibrational heat flow in a molecular lattice system consisting of a chain of particles connected to two heat baths at different temperatures, where the temperature difference between baths is oscillating in time.We derive expressions for heat currents in this system using a stochastic energetics framework and a nonequilibrium Green's function approach that is modified to treat the nonstationary average energy fluxes.We find that emergent energy storage, energy release, and thermal conductance mechanisms induced by the temperature oscillations can be controlled by varying the frequency, waveform, and amplitude of the oscillating gradient.The developed theoretical approach provides a general framework to describe how vibrational heat transmission through a molecular lattice is affected by temperature gradient oscillations.

I. INTRODUCTION
Many natural and technological systems derive their function from energy transport processes.  At e nanoscale, energy transport often involves a complex interplay and competition between multiple mechanisms and physical processes.For example, in nanotechnologies, energy transport may involve vibrational (phononic), electronic, and radiative mechanisms.5][26][27][28][29][30][31][32][33][34][35][36] Most studies of nanoscale energy transport are conducted at steady state under conditions in which the system distributions and the driving thermodynamic forces (such as a temperature gradient) are static in time.This gives rise to average energy fluxes that are constant in time.There are, however, a number of interesting and emergent energy transport properties that occur when a system is driven from steady state by a time-dependent force such as a time-varying temperature gradient.In these cases, the system will generally not reach a steady state, and will instead relax to a time-dependent nonequilibrium state.In this article, we investigate the vibrational heat flow in such a system, specifically in a model molecular lattice structure connecting two heat baths with different temperatures, where the temperature difference between baths is oscillating in time.
8][39][40][41] The core issue that induces the deviations is the interplay between diffusive and ballistic transport mechanisms, the latter being a mechanism that is not typically prominent in macroscale heat transport.
a) Electronic mail: galen.craven@gmail.comBallistic transport will generally not satisfy FL behavior, which is derived from a diffusive linear response picture.Although it is important to note that it cannot be stated generally that if a system's transport mechanism is diffusive, FL behavior will be observed.In fact, nanoscale systems that are dominated by diffusive transport often do not follow FL. 418][39] Nanoscale heat transport must therefore be treated using other theoretical principles.
The first experimental measurements of single molecule thermal conductance have recently been performed. 42,43These measurements have confirmed anomalous thermal transport behavior with respect to FL over molecular length scales and have also opened opportunities to engineer and then realize novel transport properties at the molecular level.A common experimental setup for probing single-molecule transport is a molecular junction, a system consisting of a molecular bridge connecting two electrodes.Transport through the molecular bridge, both electronic and thermal, can be induced by applying a voltage bias and/or temperature bias to the junction setup.][39][51][52][53][54][55][56][57][58] Lattice models have been instrumental in the study of heat transport and energy distribution at the nanoscale where the thermal conductance of molecules and materials can differ significantly from macroscale behaviors as discussed above.9][40]52 Typical examinations of molecular heat transport are conducted under the condition of a static temperature gradient, i.e., a temperature gradient that does not change in time.Time-dependent temperature gradient oscillations can alter a system's heat transport properties in comparison to the static case.Theoretical frameworks to understand systems with time dependent temperatures have been developed, 31,34,[59][60][61][62][63][64][65][66][67][68][69][70] and we have developed an approach to describe heat transport in such systems. 71][74][75] We have previously examined the heat conduction properties of a single particle in contact with two thermal baths, where the temperature difference between the baths is oscillating in time. 71We found that complex energy flux hysteresis and energy storage processes were induced by the temperature oscillations.2][83] Here, we extend our previous work on a single particle system by examining a molecular chain consisting of multiple interacting particles in the presence of an oscillating temperature gradient.
The remainder of the article is organized as follows: Section II contains the details of the molecular model we use to examine energy transport across an oscillating temperature gradient.In Sec.II A, the definitions and general formalism for the energy fluxes and thermal conductance in the model are presented.Section III contains derivations of the energy flux expressions for the system and the thermal baths.Section IV contains results and discussion about how the energy storage, energy release, and thermal conductance mechanisms induced by the temperature oscillations can be controlled by varying properties of the oscillating gradient.Conclusions and implications of the work are presented in Sec.V.

II. MODEL DETAILS AND FORMALISM
The model we consider is a molecular chain, i.e., a onedimensional lattice, consisting of N particles with equal masses connected by harmonic interactions.The first particle in the chain (the leftmost particle, i.e, particle 1) is coupled to a heat bath L with time-dependent oscillatory temperature T L (t) and the last particle (the rightmost particle, i.e, particle N ) is coupled to a heat bath R with time-dependent oscillatory temperature T R (t).A schematic diagram of the model is shown in Fig. 1.In order to distinguish the heat baths from the molecular chain, we will refer to the chain as the "system" and the heat baths as "baths".In the limit of constant temperatures in the baths, T L (t) = T (0) L and T R (t) = T (0) R , this model becomes a paradigmatic and well-studied model for heat transport in nanoscale systems.Here, we take the temperatures to oscillate in time, and show that these oscillations give rise to transport phenomena that are not present in the limit of a static temperature gradient.We will work in the limit that the system dynamics can be described classically.The Langevin equations of motion for this system are: where m is the particle mass (taken to be the same for every particle), x j is the displacement of particle j from its equilibrium position, k is the force constant between particles, k pin is the force constant for a pinning potential for the two edge particles, γ L and γ L are coupling constants between the system and the respective baths, and ξ L (t) and ξ R (t) are stochastic forces that obey the correlations where k B is the Boltzmann constant and T L (t) and T R (t) are the time-dependent temperatures of the left and right bath, respectively.We will consider the case where the temperatures of each bath take the specific forms: where T (0) L and T (0) R are the temperatures of the two baths in the limit of vanishing of oscillations, ∆T L and ∆T R define the amplitude of the oscillations, and ω L and ω R are oscillation frequencies.For ease of exposition, we define a temperature parameter for the system as but it is important to note that because of the oscillating gradient there is no local temperature in the chain.The friction coefficients in the model used here are static, but memory effects can be included. 57,61,62We consider here only cases in which ω L and ω R are commensurate so that the system is periodic with total period T .

A. Energy Flux Formalism
The oscillating temperature gradient in the system gives rise to time-dependent energy fluxes in the system FIG. 1. Schematic diagram of an example molecular chain with N particles connecting two heat baths with oscillating temperatures.The temperatures of the left and right baths, TL(t) and TR(t), are time-dependent and oscillate periodically with different frequencies.Because of the oscillating temperature gradient, the energy fluxes JL(t) and JR(t) in/out of the two baths are also time-dependent and oscillatory.Adjacent particles in the chain are connected by a harmonic potential with force constant k.
Using the stochastic energetics formalism of Sekimoto, 84 the energy fluxes in our model can be separated into three terms: 1. J sys is the energy flux in/out of the system 2. J L is the the energy flux associated with the left bath 3. J R is the the energy flux associated with the right bath It is important to note that in the presence of static temperature gradient under steady-state conditions, the J sys term goes to zero and the heat current through the system is defined by J L = −J R .The principal goal of this manuscript is to describe how an oscillating temperature gradient affects these energy fluxes in a harmonic lattice and to understand how those energy fluxes in turn alter thermal conductance.The energy fluxes for the system, left bath, and right bath can be defined as, where E(t) is the total energy of the system and ⟨⟩ represents an ensemble average.A more general form of energy fluxes will be developed below.
Because of the temperature oscillations, the system does not reach a steady state.Instead, the system will approach a periodic nonequilibrium state in the long-time limit (t → ∞).This means that the system energy flux J sys (t) does not go to zero when the temperatures of the baths are time-dependent.Compare this with the case in which the temperature gradient is static and the system reaches a steady state in the long-time limit.In this case, J sys (t) = 0.The implication of this is that in the periodic nonequilibrium state, the system can store and then release energy due to the temperature oscillations.These storage/release properties are prominent when the temperature oscillation frequencies are fast relative to the relaxation rates into the baths parameterized by γ L and γ R .We define the average energy storage capacity of the system through the property, where Θ is the Heaviside function and is the amount of time over an oscillation period the system heat current is positive.The metric I + sys measures the average positive system energy flux and is directly related to the amount of energy stored by the system over a period of oscillation.The I + sys and T + integrals are taken over the time interval [0, T ], i.e., over a period of total system oscillation.

III. MOLECULAR HEAT TRANSPORT
Our principal goal is to develop analytical expressions to understand heat transport through the harmonic chain in the prescence of an oscillating temperature gradient.Harmonic interactions dominate the heat transport properties in a multitude of systems, and therefore we expect that the developed theory will accurately describe a number of molecular systems.However, the inclusion of anharmonic interactions between particles can be important and may give rise to new phenomena, but is beyond scope of this work.
To derive the energy fluxes, we will work in Fourier space.Transforming the equations of motion in Eq. (1), the system of equations in frequency space becomes where x and ξ are the respective Fourier-transformed representations of x and ξ.
Nonequilibrium Green's Function (NEGF) approaches for phononic transport 85 can be employed to evaluate the energy fluxes. 58The application of the NEGF formalism to compute transport properties is a common and well-established approach to understand heat conduction in molecular structures under steady-state conditions. 2,54,57,86,87Here, because the examined system does not approach a steady state, in order to include time-dependent properties of the temperature oscillations we will need to modify the standard NEGF formalism using a bottom-up approach from the equations of motion.Developing this modified NEGF method will allow us to derive expressions for the energy fluxes.
Using matrix notation, the equations of motion for the particles in the molecular chain can be written as where XT = {x 1 , x2 , • • • , xN } is a Fourier-transformed displacement vector obtained using the relation for notational simplicity.The function G(ω) is a Green's function with conjugate adjoint G † (ω).In the Green's function, M is a diagonal mass matrix, Φ a force constant matrix extracted from Eq. ( 1), and Γ L (ω) and Γ R (ω) are spectral functions of the heat baths that connect to the molecular system.For Ohmic baths, 57 as we will consider here, the spectral functions are N × N matrices with element Γ L (ω) 00 = mωγ L and element Γ R (ω) N N = mωγ R , while all other elements in the matrices are zero.
Applying the oscillating temperatures from Eq. ( 3), the noise-noise correlations in Fourier space are where we have used the Fourier definition of the Dirac delta function and the complex exponential form for the trigonometric function to express the final forms for the correlations.
The cross correlation terms are ξL (ω) ξR (ω ′ ) = ξR (ω) ξL (ω ′ ) = 0. Defining Λ L = Γ L /ω and Λ R = Γ R /ω, the correlation functions between the stochastic noise vectors are We will now evaluate the energy fluxes for the left and right baths by combining the noise-noise correlations, the NEGF formalism, and the stochastic energetics expressions for the general energy flux expressions.Using the equations of motion in Eq. ( 11), and applying inverse Fourier transforms, we can write expressions for the correlation functions needed to evaluate the energy fluxes.In matrix notation, the left and right energy fluxes can be expressed as Therefore, deriving expressions for the energy fluxes of the baths requires evaluation of the correlation functions in Eqs. ( 16) and (17).
The first term on the RHS of Eq. ( 16) and can be evaluated as and correspondingly the first term on the RHS of Eq. ( 17) can be evaluated in a similar way yielding where we have used the identity, G T (−ω) = G † (ω).The second terms (the noise-velocity correlations) on the RHS of Eqs. ( 16) and ( 17) can be evaluated as Combining these correlation functions we now arrive at the exact expressions for the time-dependent energy fluxes of the left heat bath, the right heat bath, and the system: where Λ = Λ L +Λ R .Here, we have used the conservation of energy relation J L (t) + J R (t) + J sys (t) = 0 to deduce the system energy flux.Conservation of energy has been previously verified directly for the simpler case of a singleparticle chain. 71Equations ( 22), (23), and ( 24) are the primary analytical results of this article.They are the time-dependent energy flux expressions for the left bath, right bath, and the system in the presence of an oscillating temperature gradient.In the limit of a single particle between two oscillating heat baths, the energy fluxes can be expressed in closed form using the developed formalism (see Appendix A).

Quasistatic Limit
The quasistatic (QS) limit is the limit in which the baths oscillate slowly with respect to the system-bath couplings (ω L /γ, ω R /γ → 0, 0).In the QS limit, Eqs. ( 18) and ( 19) become These expressions can be derived by noting that terms such as (ω ± ω L ) and (ω ± ω R ) in Eqs. ( 22) and ( 23) reduce to simply ω in the QS limit, which can be shown by multiplying and dividing these terms by γ and then evaluating those terms in the QS limit.
The QS limits of the noise-velocity correlation functions in Eqs.(20) and ( 21) are where we have used the relation 25) and ( 26), we now have expressions for the heat bath energy fluxes in the QS limit: Applying the identities Γ L = ωΛ L and Γ R = ωΛ R , we can write the left energy flux as and the right energy flux as The (t) and J (QS) R (t) is the well-known expression for the static phononic transmission function.Therefore, we see that in the QS limit the heat current through the system is simply the product of a factor containing the timedependent temperature difference between baths and a frequency space integral over the same transmission function that would be obtained in the limit of a static temperature gradient.Also note that in the QS limit, (t) for all t because the system energy flux diminishes to zero, J (QS) sys = 0.This is because the energy in the system is redistributed to the baths on a much faster timescale than the temperature oscillations.

Static Limit
In the limit of vanishing temperature oscillations (∆T L , ∆T R → 0, 0), we arrive at the limit of a static (S) temperature gradient.In this limit, the temperature difference between baths is not oscillating.The energy fluxes in the static limit derived using the developed formalism are which are well-known expressions that can be obtained using Landauer's formalism for heat conduction under steady-state conditions in the classical limit.

IV. RESULTS AND DISCUSSION
To verify the accuracy of the derived theoretical expressions for the energy fluxes, we performed molecular dynamics (MD) simulations of the harmonic chain model defined in Eq. ( 1) and measured the energy fluxes in those simulations.The specific method we used is the stochastic nonequilibrium dynamical simulation approach described in Ref. 86.A comparison between the simulation results and the theoretical predictions is shown in Fig. 2. The theoretical predictions are in excellent agreement for the left, right, and system energy fluxes under the chosen system and bath parameters values.We have verified that similar agreement is observed for other system parameters and chain lengths.It is important to note that in the results shown in Fig. 2, the right bath temperature is not oscillating (i.e.∆T R = 0)-only the left bath is oscillating-but the right heat bath still participates in the heat conduction process by giving and absorbing energy from the system.This is shown in the middle panel of Fig. 2 where the right bath energy flux is oscillatory because of time-periodic energy propagation from the oscillating left bath.
The time-dependent energy fluxes for two different  chain lengths are shown in Fig. 3.The presented results provide insight into the combined effect of temperature oscillations and variation in the length of the molecular chain.Because of the oscillations, there is a timeperiodic temperature gradient ∇T (t) = ∆T (t)/N across the system.There are two key observations: (a) complex patterns arise when the temperatures of both baths are oscillating in time due to the multiple frequencies and accompanying resonances that arise in the system and (b) small differences are observed in the energy fluxes when the chain length increases from N = 2 (solid black) particles to N = 100 (dashed red) particles, suggesting that ballistic transport dominates for purely harmonic chains even in the presence of a oscillatory temperature gradient.
A natural question that arises is how heat conduction is affected by the strength of the interatomic interactions (i.e., the force constants) within the molecular chain under external temperature driving.Figure 4 illustrates the system energy flux for different force constants in a chain of N = 60 particles.The force constant in the model is denoted by k and, for ease of exposition, we use a baseline force constant value of k 0 = 250 and express changes in the force constant as multiplicative factors of this value.The force constants in the model are always the same for each bond in the molecular chain.As shown in Fig. 4, the stronger the interatomic interaction, the larger the maximum magnitude of the energy fluxes are.For example, k = 10k 0 exhibits approximately a 50% increase in energy flux magnitude compared to k = 0.1k 0 .One might naively attribute this observation to stronger correlations between the atoms in the molecule when the the force constant is larger.However, as we will demonstrate later, this interpretation does not account for the complex interplay between system-bath interactions.Figure 5 illustrates the length-dependent average energy storage capacity of the system over one period of oscillation for varying force constants.The average energy storage capacity I sys is defined in Eq. 8.The fluctuating behavior is observed for three force constants and chain lengths, reflecting the oscillatory temperatures from the baths.However, these fluctuations are not strictly in the time domain and can be nonlinear and oscillatory when other system variables are varied besides time.Several observations are of note: (a) It is not a straightforward assertion that stronger interatomic interactions result in higher conduction.For instance, specific chain lengths exhibit smaller energy storage capacities for k = 10k 0 compared to k = 0.1k 0 , although, in general, the values for k = 10k 0 are higher.(b) Smaller force constants appear to yield smaller oscillating magnitudes.This trend can be attributed to larger force constants increasing the upper limit of frequencies of hybridized modes within the chains which amplifies conduction through different channels in comparison to cases with smaller force constants.(c) There are discernible trends for all three heat current intensities to approach certain plateaus.They commence with larger fluctuations and stabilize as the chain lengths increase.This observation aligns with steady-state length-dependent heat conduction principles, where the thermal conductance of long chains becomes length-independent. 54A significant feature in this work is that the oscillations of the energy storage capacities do not vanish, but instead fluctuate about certain baseline values.
Interesting transport phenomena emerge as the driving frequencies of the baths are modified, as demonstrated in Fig. 6 which illustrates the system's average energy storage capacity over different oscillation frequencies.For simplicity, we consider the case in which only the left temperature is oscillating (i.e., ∆T R = 0).The dips in the curves signify resonant states between vibrational frequencies in the system and bath oscillation frequencies.Taking the case of k = k 0 (corresponding to k 0 /m ∼ 5ps −1 in our example), the first significant dip occurs around ω L ∼ 5ps −1 , suggesting a hybridized mode around the k 0 force constant frequency could be contributing the transport.The trends for smaller k values exhibit earlier occurrences of dips, while when k 0 is increased by a factor of ten above the baseline value, resonances under ω L = 10 values are not observed.We have confirmed that dips do occur when ω L > 10.One observation regarding higher driving frequencies is that when the bath oscillates extremely fast, the average energy storage capacities for all three force constants examined appear to converge to the same asymptotic value.This implies the energy storage capacities for the system gradually become agnostic to the force constant value.Conversely, at the opposite end of the spectrum, when the bath temperatures oscillate extremely slowly, the energy storage capacity diminishes to zero.This is a result of the vanishing system energy fluxes as the system gradually approaches the quasistatic limit.The example in Fig. 6 is presented for N = 5 particles, but similar behaviors are observed for longer chains.
Figure 7 illustrates the storage capacity as a function of the left system-bath coupling (γ L ).The right systembath coupling is held constant at γ R = 1 and the force constants inside the chain (composed of N = 5 particles) are uniformly set to k = k 0 .When the coupling on the left is extremely small, indicating that the system is nearly decoupled from the left bath, the energy storage vanishes.Conversely, for strong system-bath coupling, the capacity for the system itself to retain energy dimin-FIG.7.
Average energy storage capacity as function of system-bath coupling for a chain length of N = 5 and a force constant k = k0 = 250.The capacity is calculated using Eq. 8 and is shown in units of γkBT , a unit which is changing as γL is varied.Time is shown in units of the total oscillation period T .Parameters are m = 1, T ishes gradually because heat conduction through the system occurs quickly due to the strong coupling.The implication of this is that heat moves more quickly through the system and less energy is retained.This nonlinear pattern once again highlights the complicated interplay between all factors (dynamical and system-specific) in the transient thermal transport mechanism examined in this work, as elucidated by our analytical results.

V. CONCLUSIONS
We have examined vibrational heat transport through a harmonic molecular chain in the presence of a temperature gradient that is oscillating in time.Using a nonequilibrium Green's function approach modified to account for the time-dependent thermal gradient, we have developed a theoretical framework to describe the heat transport properties in the examined system analytically.The presented results demonstrate that the time-periodic modulation of a temperature gradient can significantly affect the heat transport and energy storage properties of nanoscale systems.
We have specifically shown that the introduction of an oscillating temperature gradient modifies the thermal conductance through a molecular lattice leading to altered heat transport, energy storage, and energy release.These modifications arise due to differences in magnitude between the system-bath relaxation rates and other frequencies in the system, for example the oscillation frequency of the temperature gradient.By tuning the properties of the oscillating gradient, such as its waveform, frequency, and amplitude, it is possible to manipulate the amount of energy stored within the system.This effect has potential applications for improved energy har-vesting and energy storage technologies, such as thermal batteries.
While the work in this article focuses on establishing a mathematical framework to describe heat transport under a time-periodic temperature gradient, the effects we have predicted could potentially be observed experimentally using current set-ups capable of measuring heat transport at the single-molecule level.][90] Investigating how anharmonicities in the lattice structure affect heat transport is a potentially important next step.Our preliminary findings suggest that when the system is asymmetric and the interaction forces between particles are anharmonic, the time-periodic modulation of a temperature gradient can give rise to emergent rectification effects.Examining these effects is a focus of our current work.Applying nonequilibrium thermodynamic relations 91 and formalisms to separate the time-dependent currents into positive and negative components [92][93][94] will provide further insight into the thermodynamic properties of the system.The present work focuses on the classical limit of temperaturemodulated heat conduction.An important future direction involves extending the model to incorporate quantum mechanical effects.This extension will be needed in order to understand the effects described in this article in regimes in which quantum effects alter transport, for example, at low temperatures.Overall, the findings in this article illustrate that nanoscale systems can be tailored to exhibit controllable heat transport properties and controlled energy storage through the temporal modulation of a temperature gradient.

FIG. 2 .
FIG. 2. Time-dependence of the energy fluxes for the left bath (top), right bath (middle), and the system (bottom) in the case of a molecular chain consisting of N = 5 particles connected to two heat baths with oscillating temperatures.In each panel, the black dashed curve is the exact analytical result and the colored line is the result generated from molecular dynamics simulations.Each energy flux is shown in units of γkBT .Time is shown in units of the total oscillation period T .Parameters are γ = 2 (γL = 1, γR = 1), m = 1, T (0) L = 1, T (0) R = 1.05, ∆TL = 0.1, ∆TR = 0, ωL = 5, ωR = 0, k = 250, kpin = 0.All parameters throughout are given in reduced units with characteristic dimensions: σ = 1 Å, τ = 1 ps, m = 10 mu, and T = 300 K.The y-axis values are running averages calculated from sets of 10 consecutive data points to reduce noise.

FIG. 4 .
FIG. 4. Time-dependence of the energy flux for a system consisting of N = 60 particles connected by harmonic forces.The different curves correspond to the different force constants shown in the figure legend.Each energy flux is shown in units of γkBT .Time is shown in units of the total oscillation period T .Parameters are γ = 1.7 (γL = 1.5, γR = 0.2), m = 1, T