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.^{1–24} At the 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. The development of theoretical methods to understand nanoscale heat transport has a rich history.^{1,2,4–6,10,25–36} Most studies of nanoscale energy transport are conducted at a 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 a 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.

At the nanoscale, heat transport properties often violate Fourier’s law (FL)—the primary physical principle that is used to understand heat transport at the macroscale.^{37–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. Ballistic 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.^{41} The physical structure of a system and its environment must often be purposefully engineered to generate FL behavior.^{37–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,43} These 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. Molecular junctions have been studied extensively and are well-used to probe molecular transport mechanisms.^{42–50}

From a theoretical standpoint, lattice models provide tractable and useful mathematical frameworks to examine molecular thermal conductance.^{4–6,37–39,51–57} 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. For example, the deviation from FL behavior at the nanoscale was predicted well before the recent experimental setups capable of measuring single molecule thermal transport were developed.^{37–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,58–69} and we have developed an approach to describe heat transport in such systems.^{70} In practical settings, temperature oscillations can be achieved using various methods, such as applying a thermal source that switches on and off at periodic intervals, using a laser to heat the system in a periodic manner, applying flash heating techniques, and/or by modulating the system pressure.^{71–74}

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.^{70} We found that complex energy flux hysteresis and energy storage processes were induced by the temperature oscillations. These effects have potential applications in, for example, thermoelectric and pyroelectric materials,^{32,75–77} thermal batteries,^{78,79} and energy transport devices.^{80–82} 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 this article is organized as follows: Sec. 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 this work are presented in Sec. V.

## II. MODEL DETAILS AND FORMALISM

The model we consider is a molecular chain, i.e., a one-dimensional 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, $TL(t)=TL(0)$ and $TR(t)=TR(0)$, 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.

*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

*γ*

_{R}are coupling constants between the system and the respective baths, and

*ξ*

_{L}(

*t*) and

*ξ*

_{R}(

*t*) are stochastic forces that obey the correlations

*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 following specific forms:

*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

^{56,60,61}We 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

^{83}the energy fluxes in our model can be separated into three terms as follows:

*J*_{sys}is the energy flux in/out of the system.*J*_{L}is the energy flux associated with the left bath.*J*_{R}is the energy flux associated with the right bath.

*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 article 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 follows:

*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.

*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 following property:

## III. MOLECULAR HEAT TRANSPORT

Our principal goal is to develop analytical expressions to understand heat transport through the harmonic chain in the presence 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 the scope of this work.

Harmonic interactions govern an important and fundamental characteristic of molecular heat transport—ballistic transport, a phenomenon in which thermal energy propagation through chains or lattices is not influenced by scattering effects. Ballistic transport leads to thermal conductance values that remain constant over increasing system lengths,^{5,56,84} a deviation from Fourier’s law behavior.^{38–41} Recent computational and experimental investigations have confirmed this behavior in various molecular systems, such as alkane-based chains and benzene-based ring molecules.^{22,42,43,85,86} Anharmonic interactions have been shown to alter heat conduction in certain molecular systems resulting in deviations from ballistic behavior.^{9} This generates transport effects (e.g., thermal rectification^{19,30,41,87–92}) that are not present in harmonic systems.

*x*and

*ξ*.

Nonequilibrium Green’s Function (NEGF) approaches for phononic transport^{93} can be employed to evaluate the energy fluxes.^{57} The 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,9,54,56,86,94} Here, 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.

**Ξ**(

*ω*) =

**Ξ**

_{L}(

*ω*) +

**Ξ**

_{R}(

*ω*) for notational simplicity. The function

**(**

*G**ω*) is a Green’s function with conjugate adjoint

*G*^{†}(

*ω*). In Green’s function,

**M**is a diagonal mass matrix,

**Φ**is 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,

^{56}as we will consider here, the spectral functions are

*N*×

*N*matrices with element

**Γ**

_{L}(

*ω*)

_{00}=

*mωγ*

_{L}and element

**Γ**

_{R}(

*ω*)

_{NN}=

*mωγ*

_{R}, while all other elements in the matrices are zero.

**Λ**

_{L}=

**Γ**

_{L}/

*ω*and

**Λ**

_{R}=

**Γ**

_{R}/

*ω*, the correlation functions between the stochastic noise vectors are

**G**

^{T}(−

*ω*) =

**G**

^{†}(

*ω*). The second terms (the noise–velocity correlations) on the RHS of Eqs. (16) and (17) can be evaluated as

**Λ**=

**Λ**

_{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 single-particle chain.

^{70}Equations (22)–(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 the closed form using the developed formalism (see Appendix A).

### A. Quasistatic limit

*ω*

_{L}/

*γ*,

*ω*

_{R}/

*γ*→ 0, 0). In the QS limit, Eqs. (18) and (19) become

*ω*±

*ω*

_{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.

**G**

^{†}−

**G**= 2

*i*

**G**

^{†}

**ΓG**.

**Γ**

_{L}=

*ω*

**Λ**

_{L}and

**Γ**

_{R}=

*ω*

**Λ**

_{R}, we can write the left energy flux as

**G**

^{†}(

*ω*)

**Γ**

_{L}

**G**(

*ω*)

**Γ**

_{R}] integrand in $JL(QS)(t)$ and $JR(QS)(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 time-dependent 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. In addition, note that in the QS limit, $JL(QS)(t)=\u2212JR(QS)(t)$ for all

*t*because the system energy flux diminishes to zero, $Jsys(QS)=0$. This is because the energy in the system is redistributed to the baths on a much faster timescale than the temperature oscillations.

### B. Static limit

*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

^{2,56,86}

## 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 insights into the combined effect of temperature oscillations and variation in the length of the molecular chain. Because of the oscillations, there is a time-periodic 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 (black solid) particles to *N* = 100 (red dashed) particles, suggesting that ballistic transport dominates for purely harmonic chains even in the presence of an 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* = 10*k*_{0} exhibits approximately a 50% increase in energy flux magnitude compared to *k* = 0.1*k*_{0}. One might naively attribute this observation to stronger correlations between the atoms in the molecule when 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* = 10*k*_{0} compared to *k* = 0.1*k*_{0}, although, in general, the values for *k* = 10*k*_{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.^{9} A 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* ∼ 5 ps^{−1} in our example), the first significant dip occurs around *ω*_{L} ∼ 5 ps^{−1}, suggesting that a hybridized mode around the *k*_{0} force constant frequency could be contributing to 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 that 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 system–bath 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 diminishes 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 harvesting 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 setups capable of measuring heat transport at the single-molecule level. Specific systems that may exhibit the phenomena we have predicted include alkanedithiol chains and polymeric molecular systems with metallic substrates in molecular junction setups under periodic temperature oscillations.^{42,92,95,96}

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^{97} and formalisms to separate the time-dependent currents into positive and negative components^{98–100} will provide further insights into the thermodynamic properties of the system. The present work focuses on the classical limit of temperature-modulated 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.

## ACKNOWLEDGMENTS

We acknowledge support from the Los Alamos National Laboratory (LANL) Directed Research and Development (LDRD) funds. This research was performed in part at the Center for Nonlinear Studies (CNLS) at LANL. The computing resources used to perform this research were provided by the LANL Institutional Computing Program.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Renai Chen**: Investigation (equal). **Tammie Gibson**: Investigation (equal). **Galen T. Craven**: Investigation (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: HEAT CURRENT DERIVATION FOR A SINGLE-PARTICLE (*N* = 1) CHAIN

The developed theoretical framework can be applied to generate closed-form expressions for the energy fluxes in the case of a single particle bridging two heat baths with oscillating temperatures. These results can then be compared and validated using results obtained previously through an approach that uses direct calculation of the time correlation functions in the energy flux expressions.^{70}

*γ*=

*γ*

_{L}+

*γ*

_{R}and effective noise $\xi \u0303(\omega )=\xi \u0303L(\omega )+\xi \u0303R(\omega )$ terms. Note that the noise–noise correlation functions are the same as in Eqs. (13) and (14). The equation of motion can be written as

*G*in the following derivation [compare with Eq. (11) where the mass is moved to the RHS].

*J*

_{L}(

*t*) +

*J*

_{R}(

*t*) +

*J*

_{sys}(

*t*) = 0 to obtain

*J*

_{sys}. The energy flux expressions derived here using the NEGF approach are in agreement with the results we have derived previously by directly evaluating the correlation functions for a single Brownian particle in Ref. 70.

### APPENDIX B: DERIVATION FOR THE RELATION USED IN EQ. (A15)

*G*(

*ω*), we have

*G*−

*G*

^{−}= −2

*iωγ*|

*G*(

*ω*)|

^{2}, we have arrived at Eq. (A15),

## REFERENCES

*Mathematical Physics*