Fluctuating charge distribution of one molecule can interact with the fluctuating charge distribution of another proximate molecule so as to give rise to more energetically advantageous fluctuations. This interaction is responsible for fundamentally important phenomena such as London’s dispersion forces, and that relevant for this work–near-field heat transfer (NFHT). In order to model the dominant dipolar terms of these fluctuating charge interactions, and thus the ensuing phenomena, it is possible to use a Drude oscillator model where the electric charge distribution in the atom or molecule undergoes an oscillatory displacement in response to an oscillatory electric field. This has been successfully employed to model among others the London dispersion forces between molecules [1], van der Waals interaction between atom and a surface [2] and between two surfaces [3].

For this study we analyze NFHT between two closely spaced molecules. The study of NFHT has gained prominence in recent years due to theoretical predictions [4–6] and subsequent experimental confirmations [7–14] of heat transfer between vacuum separated objects exceeding Planck’s blackbody predictions by several orders of magnitude. This has led to development of multitude of novel applications such as near-field thermophotovoltaics [15, 16], thermal diodes [17], transistors [18, 19], rectifiers [20] and modulators [21]. A more exhaustive list of applications can be found in any recent review on this topic such as in Refs. [22–24]. Experimental advancements has lead to measurement of NFHT with conductances of the order of 200 pWK⁻¹ for configurations such as that between STM tip over substrate [14, 25]. Recently it has been shown that thermal conductance in a single molecule junction between two heaters maintained at different temperatures can be measured with resolution as low as 2 pWK⁻¹ [26]. A similar configuration can be adopted to measure thermal conductance between two molecules (or two nanoparticles in general) but now with vacuum gap between them. As the molecules are brought closer, nonradiative dipole coupling will result in a gap-dependent heat flux which we attempt to quantify in this theoretical study.

In our study we adopt coupled Drude oscillators to model the interaction between the molecules and estimate NFHT between them. The advantage of using this model stems not only from the fact that this picture offers mathematical simplicity and hence the physics is transparent but also from that it establishes a common theoretical framework to analyze the two phenomena of dispersion forces (for which the theory has been well developed) and NFHT both of which have common origin in the interaction between fluctuating charge distributions. This study also enables us to predict the effect of orientation factor of two dipoles on near-field heat transfer which has so far not been taken into account in the fluctuational electrodynamics framework traditionally used to analyze near-field heat transfer. Another added advantage is that the results from Drude oscillator model can be easily extended to estimate interaction between molecules by including either the tabulated oscillator strengths [27] for the respective molecules, or by calculating them from their relations with experimental measurements of fluorescence and absorption spectra [28]. Recently [29], have proposed an alternative method to compute the heat transfer between molecules which requires computation of Green’s function for an approximate geometry of the molecule via boundary element or finite-difference methods, and the molecular susceptibilities obtained from density functional theory. In contrast, the procedure outlined in this manuscript proposes an alternative approach: to make use of experimentally determined fluorescent and emission spectra of the molecules (which inherently depend on the shape, size and composition of the molecule) to arrive at an expression for the near-field heat transfer.

The procedure for detailing the near-field heat transfer between two molecules is as follows: we first derive the eigenmodes for the case of two interacting dipoles in Dipole-Dipole Coupling for Near-Field Heat Transfer and by comparing these with the corresponding forms in the coupled harmonic oscillator model shown in Coupled Harmonic Oscillator Model we get the expressions for the parameters that can be substituted in the expression for heat transfer mediated by a coupled harmonic oscillator system between two heat baths maintained at two different temperatures. We show that results obtained from such a trivial substitution conforms to the predictions from fluctuational electrodynamics (FE) theory. The NFHT between two dipoles is then extended to predict that between two molecules by incorporating their oscillator strengths. This procedure has been previously applied to estimate the NHFT between two nanoparticles [30] and two planar surfaces [31, 32]. The model has the added advantage that it can be extended to cases where the local thermal equilibrium condition considered in FE is not valid - such as estimating dynamic heat transfer between objects [30, 33]. The procedure detailed in Dipole-Dipole Coupling for Near-Field Heat Transfer and Coupled Harmonic Oscillator Model is similar to that detailed in Ref. [32] and has been included for the purpose of completion and providing additional details in the relevant calculations. In addition the results shown in Ref. [32] differed from FE predictions by a constant factor. In this work we account for this discrepancy by taking into account degeneracy in eigenmodes and in NFHT Between Molecules we discuss the effect of orientation factor between dipoles on NFHT between molecules.

Consider two transient dipoles d i i = 1,2 which are self-sustaining due to interactions from each other’s electric fields: d 1 = α 1 ( ω ) E 2 ; d 2 = α 2 ( ω ) E 1 (1) where α i ( ω ) denotes the polarizability of the corresponding dipoles. The electric field emitted by a dipole contains [34] both components which propagate over distances longer than the wavelength of radiation (termed as “far-field” and which varies as 1 / | r | 2 and 1 / | r | where r is the displacement vector from the position of the dipole) and components which decay exponentially over distances smaller than the wavelength of the radiation (termed as “near-field” and which varies as 1 / | r | 3 ). Since we are interested only in the near-field effects, taking only the terms of the electric field from a dipole which varies as 1 / | r | 3 we have: E 1,2 = Re [ ( 3 r ^ r ^ . d 1,2 | r | 3 − d 1,2 | r | 3 ) e i k . r 4 π ϵ ] (2)

Here r ^ is the unit vector in the direction of r and k is the wavevector of emitted radiation, and ϵ is the dielectric permittivity of the medium where the dipoles are located. In the limit of small gaps between the dipoles ( k . r ≪ 1 ) we have: e i k . r ≈ 1 . The expressions for the dipoles in Eq. 2 can then be written as: d 1 = α 1 ( ω ) R ↔ d 2 ; d 2 = α 2 ( ω ) R ↔ d 1 ; (3) where, if we consider without loss of generality that the two dipoles are separated along z-axis, R ↔ = 1 4 π ε r 3 ( − 1 0 0 0 − 1 0 0 0 + 2 ) (4)

Eliminating one of the dipoles in Eq. 3 we find two conditions for non-zero solutions to exist α 2 ( ω ) α 1 ( ω ) ( 1 4 π ϵ r 3 ) 2 = 1 and α 2 ( ω ) α 1 ( ω ) ( 2 4 π ϵ r 3 ) 2 = 1 . If the two dipoles have polarizability given by the standard Lorentzian form [34]: α 1,2 ( ω ) = e 2 / m ω D 2 − ω 2 − i ω δ (5) with ω D being the natural frequency of the Lorentz oscillator and δ denoting the spontaneous decaying rate, we obtain the expressions of the eigenmodes from the two conditions as: ω ± L = ω D 1 ± e 2 2 π ϵ d 3 m ω D 2 − i δ 2 ; ( with degeneracy  1 ) (6) and ω ± T = ω D 1 ± e 2 4 π ϵ d 3 m ω D 2 − i δ 2       ( with degeneracy      2 )   (7) where, the superscripts L and T, standing for longitudinal and transverse modes, have been included to differentiate between the two possible solutions with different degeneracies. Physically these two solutions denote the possibility of the dipoles being in either “head-to-tail” configuration or “side-by-side” configuration respectively [35], also termed “J-aggregate” and “H-aggregate”, respectively [36]. The Lorentz form of polarizability has been assumed in Eq. 5 since it allows us to readily extend the analysis to estimate the NFHT between two molecules by including the appropriate spectral oscillator strength as shown in NFHT Between Molecules. By replacing this expression of polarizability with other forms applicable for larger objects, such as nanoparticles, one can trivially extend this analysis to estimate heat transfer between such objects. It should also be possible to generalize this method when there are several interacting dipolar structures such as that seen in hybrid graphene nanostructures [37].

The equations of motion of two coupled harmonic oscillators of unit mass, with same natural frequency ω 0 , damping constant γ, and coupling constant ω 0 2 ξ , with one of the oscillators connected to a heat bath at temperature T and the other to a heat bath at temperature 0 K, as shown in Figure 1, is given by [31–33, 38]: x 1 ″ ( t ) + ω 0 2 x 1 ( t ) + γ x 1 ′ ( t ) + ω 0 2 ξ x 2 ( t ) = f ( t , T ) x 2 ″ ( t ) + ω 0 2 x 2 ( t ) + γ x 2 ′ ( t ) + ω 0 2 ξ x 1 ( t ) = 0 (8) where, the forcing function f ( t , T ) is the rapidly fluctuating Langevin force generated by the heat bath.

The equations of motion in frequency space reduce to: ω 2 X 1 = ω 0 2 X 1 − i γ ω 0 X 1 + ω 0 2 ξ X 2 + F ( ω , T ) ; ω 2 X 2 = ω 0 2 X 2 − i γ ω 0 X 2 + ω 0 2 ξ X 1 (9) where X 1,2 and F ( ω , T ) are the Fourier coefficients for x 1,2 and f ( ω , T ) respectively. The natural frequencies for this system of equations are: ω 1,2 = ω 0 ( 1 ± ξ ) − i γ 2 (10) where it is assumed that γ / ω 0 ≪ 1 . For the forced oscillations [for finite F ( ω , T ) we get the solution matrix X ( ω ) to be of the form: X ( ω ) = 1 2 ( 1 ω 2 − ω 2 2 + 1 ω 2 − ω 1 2 1 ω 2 − ω 2 2 − 1 ω 2 − ω 1 2 ) F ( ω , T ) (11) where expressions for the natural frequencies ω 1 and ω 2 are given in Eq. 10. The spectral density of the forcing function, S F ( ω , T ) , which is proportional to | F ( ω , T ) | 2 can be found by equating the thermal energy content of the oscillator in the absence of coupling with that predicted from Bose-Einstein distribution as detailed in Ref. [32]. We find this to be: S F ( ω , T ) = γ π ℏ ω ( e ℏ ω k B T − 1 ) (12)

The net heat transferred in this coupled harmonic oscillator (CHO) system can be found by observing that in steady state the rate of heat transferred to the second oscillator via coupling with the first oscillator is equal to the decay in the second oscillator via damping: P CHO = − ω 0 2 d dt 〈 x 2 2 ( t ) 〉 = γ ω 0 2 〈 x 2 2 ( t ) 〉 (13)

Since we can relate 〈 x 2 2 ( t ) 〉 to | X 2 ( ω ) | 2 using power spectral density S F ( ω , T ) as: 〈 x 2 2 ( t ) 〉 = ∫ 0 ∞ | X 2 ( ω ) | 2 S F ( ω , T )   d ω , we obtain from Eqs. 11–13 noting that the integrand in 〈 x 2 2 ( t ) 〉 is sharply peaked around ω ≈ ω 0 : P CHO ≈ ℏ ω 0 e ℏ ω 0 / ( k B T ) − 1 ξ 2 ω 0 2 γ / 2 γ 2 + ξ 2 ω 0 2 (14)

We now ask: what are the equivalent expressions for the parameters ω 0 , ξ, and γ suitable for modeling the dipole-dipole interaction as given by Eq. 1? This can be found by comparing the forms of eigenmodes in Eqs. 6, 7, 10 from which we get: ω 0 → ω D ;     ξ L → e 2 2 π ϵ d 3 m ω D 2 ;     ξ T → e 2 4 π ϵ d 3 m ω D 2 ;     γ → δ (15)

Substituting the coupling parameters from above we get: P L = Θ ( ω D , T ) ( δ e 4 e 4 + 4 π 2 ϵ 2 δ 2 d 6 m 2 ω D 2 ) ; (16) P T = Θ ( ω D , T ) ( δ e 4 e 4 + 16 π 2 ϵ 2 δ 2 d 6 m 2 ω D 2 ) (17) where, Θ ( ω D , T ) = ( ℏ ω D ) / [ e ℏ ω D / ( k B T ) − 1 ] . For the frequencies that we are interested in i.e., ω D ≈ 10 15 rad s − 1 and gaps of the order 10⁻⁹ m we will have e 4 ≪ 4 π 2 d 6 δ 2 m 2 ω D 2 ϵ 2 , so that we get the expression for heat flux P between two dipoles is: P = P L + 2 P T = Θ ( ω D , T ) 3 e 4 16 π 2 d 6 δ m 2 ω D 2 ϵ 2 (18)

We note here in passing that for ω D , d, δ parameters such that e 4 ≫ 4 π 2 d 6 δ 2 m 2 ω D 2 ϵ 2 Eq. 16 predicts that the power transferred between the two dipoles will be independent of the gap between them. However, at such small gaps, other effects not considered in this work such as nonlocal nature of polarizability will also have to be taken into consideration.

Here we use the expression of heat flux between two dipoles derived in Eq. 14 to predict the heat transfer between two molecules for the general case when resonant frequencies of the two molecules are different. For the case of two oscillators with different resonant frequencies ω D and ω A the equivalent expression for the eigenmodes in Eq. 10 is given by: ω 1,2 = ( ω D 2 + ω A 2 ) / 2 1 ± 1 − ( ω A ω D ( ω D 2 + ω A 2 ) / 2 ) 2 ( 1 − ξ 2 ) + i γ ( ω A + ω D ) / 2 2 ( ω A 2 + ω D 2 ) / 2 (21) which, for cases where ω D ≈ ω A we will have: ω 1,2 ≈ ω D + ω A 2 1 ± ξ + i γ 2 (22)

Thus, as in Eq. 15, the equivalent expression for heat flux between molecules with natural frequencies ω D ≈ ω A can be found by replacing ω 0 → ( ω D + ω A ) / 2 . Employing this in Eq. 14, and noting that due to the probability distribution of the electron cloud we will have a spectrum of equilibrium positions [41] we can now estimate the near-field heat transfer between two molecules by including the respective spectral oscillator strength F 1 ( ω D ) and F 2 ( ω A ) as: P = 3 e 4 16 π 2 ϵ 2 d 6 δ m 2 ∫ 0 ∞   ∫ 0 ∞ Θ ( ω D , T ) F 1 ( ω D ) F 2 ( ω A ) d ω D d ω A ( ( ω D + ω A ) / 2 ) 2 (23)

Values of F 1 ( ω D ) and F 2 ( ω A ) for different molecules can be found from tabulated data [27] or can be calculated from their relations to the experimentally determined absorption cross section σ A ( ω ) of the acceptor molecule, the normalized fluorescent spectrum f D ( ω ) , and the quantum yield Φ D of the donor molecule as [28, 42]: F 1 ( ω ) = 3 c 3 ( 4 π ϵ 0 ) Φ D m e f D ( ω ) 2 e 2 n ω 2 τ D (24) and F 2 ( ω ) = 3 c n ( 4 π ϵ 0 ) m e σ A ( ω ) π e 2 (25)

Here, c is the velocity of light in vacuum, ε 0 is the permittivity of free space, m the mass of electron, and τ D the fluorescence life time, and n is the refractive index of the surrounding medium.

We now consider the case of heat transfer between molecules when the relative orientation of dipole moments is fixed (possibly due to the molecules not being free to rotate) and indicated by the orientation factor κ. This situation is frequently encountered in the analysis of Förster resonance energy transfer (FRET) between two adjacently placed molecules. In fact, Förster in his original treatize on FRET [42] analyzed the interaction between molecules using two drude oscillators (an english version of the derivation can be found in Ref. [28]). We now consider the effect of this orientation factor on NFHT between the molecules.

Since only the component of electric field of one dipole along the other dipole is responsible for work done and hence the heat flow, we consider only the component of the electric field given in Eq. 2 along the driven dipole [43–45]: E 1,2 = d 1,2   κ 4 π ϵ d 3 (26) where κ = cos θ 12 − 3 cos θ 1 r cos θ 2 r with θ 12 being the angle between the dipoles d 1 and d 2 , and θ i r (i = 1,2) being the angles between d i and r ^ . Here, ϵ m is the relative permittivity of the medium separating the two molecules. Eliminating one of the dipoles in Eq. 1 as done previously in Dipole-Dipole Coupling for Near-Field Heat Transfer gives us the condition for nonzero solutions to exist as: 1 − α 1 ( ω ) α 2 ( ω ) ( κ 4 π ϵ d 3 ) 2 = 0 (27)

For two similar molecules of polarizability α ( ω ) of the form given in Eq. 5 we get the two solutions by solving α ( ω ) κ / ( 4 π ϵ d 3 ) = ± 1 which gives us: ω ± = ω D 1 ± e 2 κ 4 π ϵ d 3 m ω D 2 + i δ 2 (28)

Comparison with Eq. 10 gives ω 0 → ω D ;     ξ → e 2 κ 4 π ϵ d 3 m ω D 2 ;     γ → δ (29)

Using this in the expression for heat flux between two dipoles in Eq. 18 gives: P = Θ ( ω D , T ) e 4 κ 2 32 π 2 ϵ 2 δ d 6 m e 2 ω 0 2 (30)

The equivalent expression for heat flux between two molecules taking into account the dipole orientation factor κ in Eq. 23 will then be: P = e 4 κ 2 32 π 2 ϵ 2 d 6 δ m 2 ∫ 0 ∞   ∫ 0 ∞ Θ ( ω D , T ) F 1 ( ω D ) F 2 ( ω A )     d ω D   d ω A ( ( ω D + ω A ) / 2 ) 2 (31) where expressions for F 1 ( ω D ) and F 2 ( ω D ) are given in Eqs. 24, 25, respectively. For the case when all dipole orientations are equally likely, we have to add up the contributions of heat flux from all the channels: two from the transverse mode ( κ = 2 ) and one from the longitudinal mode ( κ = − 1 ) and the expression for heat flux reduces to that between two isotropic molecules given in Eq. 18. Recently it has been shown that orientation factor can be engineered to control FRET [46]. Engineering such a configuration while measuring NFHT will confirm the role of the orientation factor in the heat transfer.

For demonstration, we now calculate the thermal conductance, defined as G = d P d T | T = 300 where, p is given by Eq. 30 for heat flux between two dipoles, and given by Eq. 31 for heat flux between two molecules. Consider first the case of two dipoles whose polarizability is given by Eq. 5 separated by a gap of 2 nm in vacuum and oscillating with natural frequency ω D = 1.95 × 10 14 rad/s which corresponds to the peak of the black body spectrum at 300 K. Considering the alignment of the dipoles to be constrained such that κ = 2 , and a value of line width δ = 0.025 ω D , we obtain value of conductance between two dipoles G dip = 3.9   × 10 − 8 WK⁻¹. This value of conductance between two dipoles at T = 300 K is consistent with that from FE predictions [47]. For calculating thermal conductance between two molecules we have to use the expression for p given by Eq. 31 for which the experimental parameters σ A ( ω ) , f D ( ω ) , Φ D , τ D , has to be known. These can be obtained from catalog of spectra of molecules such as that in Ref. [48]. For this work we consider the general case where the shape of the emission and absorption lines can be estimated to be Lorentzian i.e., we take f D ( ω ) = δ / π ( ω D − ω ) 2 + δ 2 , σ ( ω A ) = δ / π ( ω D − ω ) 2 + δ 2 π e 2 2 ε 0 m c . The quantum yield Φ D and the fluorescent lifetime τ D are taken to be 0.5 and 1.5 ns, respectively. The exact values for any specific molecule is not expected to vary significantly from these estimates. Retaining the values of ω D = ω A = 1.95 × 10 14 rad/s, δ = 0.025 ω D , κ = 2 , n = 1 , and d = 2 nm we obtain the conductance between the two molecules as: G mol = 4.4 × 10 − 5 WK⁻¹, a value which is well within the limits of current experimental capabilities as detailed in Introduction.

We have shown here the framework that can be used for analyzing the NFHT between two molecules. This method employs the expression of power flow between two coupled Drude oscillators, and the necessary coupling parameters required to model the near-field dipole-dipole interactions. The advantage of adopting this approach is that the expression for power flow between the two Drude oscillators can be linked trivially to the NFHT between the two molecules via the experimentally determinable absorption and fluorescent spectra of the molecules as shown in NFHT Between Molecules. The effect of orientation factor on the NFHT between molecules is also considered. We have confirmed that this model in the dipole limit tallies with the results of NFHT between two dipoles as predicted from fluctuational electrodynamics principles and the discrepancy of a constant factor which was reported in an earlier work [32] has been accounted for by considering degeneracy of eigenmodes. It is expected that the expression of NFHT between molecules derived in this work will be used for comparison with experimental measurements. Further, by modifying the coupling term in Eq. 15 this method can be extended to include non-local effects arising from nuclear motion [29], an exercise which will be carried out in a follow-up article.