The theory of photon-mediated heat transfer between macroscopic objects in close proximity to each other and separated by a vacuum gap has traditionally been treated using the macroscopic Rytov’s fluctuational electrodynamics theory, which assumes local thermodynamic equilibrium in the bodies in question , , , . This heat transfer comprises contributions from both long-range radiative modes and near-field evanescent and surface modes , . When thermally excited, contributions from surface modes, which are electromagnetic eigenmodes of the surface and are characterized by their field decaying exponentially on either side of the interface, dominate the heat transfer between the surfaces when the gap is less than the thermal wavelength of radiation. This is primarily due to a peak in the density of electromagnetic states at such frequencies where these surface modes are resonantly activated as evidenced from the dispersion relation for these modes . In particular, for this effect to be prominent at room temperature, the surfaces should be made up of a polar crystal such as silicon carbide or silica, which supports surface-phonon polariton modes in the infrared wavelength of about 10 μm and can thus be thermally excited at these temperatures.
In general, the resonant excitation of surface modes plays an important role in several phenomena and applications, including decreased lifetime of molecules close to metal surfaces , surface-enhanced Raman spectroscopy , thermal near-field spectroscopy , , concept of perfect lens , and thermal rectification . The study of coupling of surface modes across surfaces is significant, as it plays an important role not only in heat transfer but also in the van der Waals and Casimir forces between them , . A coupled harmonic oscillator (HO) description for the heat transfer between nanoparticles due to the coupling of surface modes was arrived at by Biehs and Agarwal , and estimates for both dynamic and steady-state heat transfer values were arrived at. Barton  has considered the heat flow between two HOs using Langevin dynamics and has extended this model to planar surfaces but has limited his description to the steady-state heat flow. Yu et al.  have recently analyzed the dynamics of radiative heat exchange between graphene nanostructures using fluctuational electrodynamics principles and have observed thermalization within femtosecond timescales. This ultrafast heat exchange due to the time-varying temperatures of the nanostructures is a result of the low heat capacity of the graphene nanostructures and the coupling of the large plasmonic fields. A similar analysis for the radiative heat transfer between graphene and a hyperbolic material has been carried out by Principi et al. , where they observed thermalization in picosecond timescales. In this paper, we modeled the dynamic heat transfer contribution from the coupling of surface modes across two dielectric planar surfaces using the master equation description of two coupled HOs interacting with their respective heat baths and compared the steady-state results to those obtained from fluctuational electrodynamics principles available in the literature , . Our work differs from that of Refs. ,  in that we are interested in analyzing how the coupling between the surface modes and the resultant heat transfer between the two surfaces, which are maintained at fixed temperatures, relaxes to steady state.
The paper is arranged as follows. In Section 2, the results of the dynamics of heat transfer between two coupled HOs, each in contact with a heat bath, are summarized. Although the theory has been described in detail in Ref. , for the sake of completion, the main results have been reproduced here with added details. Such a system has also been analyzed previously in the context of analyzing the dependence of mean interaction energy on the temperatures of heat baths  and entanglement between two particles . In Section 3, the theory of dynamics of heat transfer between two coupled HOs is extended to that between two half-spaces by analyzing the coupling between two interacting planar surface modes using Maxwell’s equations. In Section 4, the numerical values for the heat flux derived in Sections 2 and 3 are plotted for the particular case of two silicon carbide half-spaces and compared to calculations from fluctuational electrodynamics principles.
2 Heat transfer dynamics between two HOs
A surface polariton located at the interface between a half-space located at z>0 and vacuum will have field of the form exp(ik.r – iωt+ikzz), where the in-plane component k=(kx, ky) and the z-component kz of the wavevector are related as where c is the velocity of light, ω is the frequency of the planar wave in the vacuum gap of this system, and k is limited by |k|>ω/c. The surface polariton exists for the (k, ω) pair that satisfies the well-known dispersion relation  so that the frequency ω0 of the surface polariton is characterized by wavevector k. Thus, one can characterize the surface excitations in terms of oscillators with frequency ω0(k), complex amplitude a(k), and half linewidth γ(k), with both ω0(k) and γ(k) determined from the dispersion relation.
When two half-spaces are far apart, then each half-space consists of oscillators with frequencies ω0(k), which are in thermal equilibrium at the temperature of the half-space. When they are brought close to each other as shown in Figure 1, then the surface polariton modes of half-space A interact with the surface polariton modes of half-space B. However, it should be borne in mind that wavevector k is conserved for fields across planar interfaces. This implies an effective coupling of the form:(1)
where b(k) is the complex amplitude of the surface polariton on half-space B, g(k) is the coupling constant (rad s−1), and we neglect nonresonant contributions to the coupling (i.e. we use the rotating wave approximation). Thus, mode a(k) couples to mode b(k) only. There are no coupling terms of the form a†(k1)b(k2), k1≠k2. Thus, we can consider coupling between two oscillators a(k) and b(k) for a fixed k and at the end of calculation sum over all modes. Note also that in the nonretarded limit the surface modes occur for all values of k. This coupling results in a heat flow due to the interaction between the modes, which is a dynamical process and can be described using the master equation approach. The structure of the weakly coupled master equation for two interacting oscillators [with amplitudes a(k) and b(k)] is well known  and is given by(2)
Here, ρS is the density matrix for the two oscillator system. The operators a, a†, b, and b† satisfy Boson algebra with the symbol [..] denoting the commutation operator, and for brevity, we drop the argument k. We have assumed that the two half-spaces have identical dielectric properties. Half-space A is at temperature T so that n1 is the excitation of mode a, i.e. n1=1/[exp(ħω0/(kBT)–1)], where kB is the Boltzmann constant. Half-space B is at zero temperature. We also take g<<ω0, which is later shown to be a valid assumption for the system under consideration. Higher-order terms in the master equation would be necessary when this assumption is not satisfied . This master equation approach is valid only for oscillators where γ/ω0<<1 and can be used to describe dynamics for time scales t>>1/ω0.
From this, and the relating expectation value of an operator 〈G〉 to the density matrix ρS using the standard relation ∂〈G〉/∂t=Tr(∂ρS/∂tG), we obtain the rate of change of the excitation of surface polariton in B to be:(3)
The first term on the right-hand side represents the change in the excitation of surface polariton in B due to the flow of thermal excitation from the surface polariton in A to that in B due to the temperature gradient. This enables us to identify the heat transferred (W) to surface polariton in B as(4)
To find 〈P〉(t), we follow the procedure used to obtain Eq. (3) to get the dynamic matrix equation:(5)
with x=(〈a†a〉, 〈b†b〉, i(〈b†a〉−〈a†b〉))T, a=(2γn1, 0, 0)T and
The value for i(〈b†a〉−〈a†b〉) in Eq. (4) from which all time-dependent and steady-state characteristics of the heat transfer can be obtained is got by solving the dynamic matrix equation in Eq. (5) with the initial conditions: 〈a†a〉|t=0=n1; 〈b†b〉|t=0=0; i〈b†a−a†b〉|t=0=0, i.e. we consider the interaction to switch on at time t=0 and analyze the dynamical evolution of the system. We thus get the expression for heat transfer from Eq. (4) as(6)
3 Heat transfer between two half-spaces
In this section, we discuss how the different parameters in the master equation can be obtained from the nature of surface polaritons for a system of two dielectric half-spaces separated by a vacuum gap and subsequently calculate the heat transfer for this system. Consider first a system of two coupled oscillators with natural frequency ω0, half linewidth γ, coupling constant and displacements x1(ω), x2(ω), where which are related as(7)
The eigenfrequencies for such a system are given by(8)
To find the equivalent coupling constant ξ between surface modes, we find the eigenmodes of the configuration of two flat surfaces separated by a vacuum gap using Maxwell’s equations and compare the expression to that of Eq. (8) for two coupled HOs. Consider two half-spaces separated by a gap of width d occupying the regions z<−d/2 and z>d/2 as shown in Figure 1. From the dispersion relation for surface polaritons ω(k) and close to the surface polariton frequency, the surface mode is seen to acquire an electrostatic character with |k|→∞ . In this electrostatic limit, satisfying the continuity of perpendicular component of displacement field gives the condition for surface modes as(9)
where ε(ω) is the dielectric function of the half-space. Taking the Lorentz model for the dielectric function:(10)
and solving for ω in Eq. (9) in the low-loss limit Γ → 0, and with some algebra, we obtain the eigenfrequencies of the coupled surface modes to be of the form:(11)
In the limit of large gaps |k|d→∞, where ω0 is the surface polariton frequency for a single half-space given by and the coupling parameter ξ′(k, d) ≈ 0. For the case when ε∞=1, these expressions reduce to the simple forms: and . This method of finding the eigenmodes of the surface modes in the electrostatic limit is well known and has been previously employed, for example, in the context of deriving the van der Waals force  and the steady-state van der Waals heat transfer . By comparing Eqs. (11) and (8), we can relate the parameters of the HO model with those of the surface modes:(14)
We can thus use the results from Section 2 to find the heat transfer between the two flat surfaces due to the coupling between the surface modes. As noted earlier, we add the contributions from all the HO modes (labeled by k), observing that heat transfer does not result in change in momentum and that only surface modes of the same in-plane wave-vector component interact across the two half-spaces. The heat flux between two half-spaces P(t, d) (W m−2) is then given by(15)
where 〈P〉(t, k, d) is the rate of heat transfer (W) between two HOs obtained from Eq. (4). The relation between g in Eq. (4) and ξ in Eq. (8) can be derived by comparing the interaction Hamiltonians in the quantum mechanical picture and the classical model. In the classical picture, the potential energy of the Hamiltonian of the coupled HO system is given by
where the interaction term is given from The equivalent quantum mechanical model can be obtained by substituting the displacements x1 and x2 in terms of the commutator operators a, a†, b, and b† using the standard relations and Including only terms resulting in photon exchange, we get
Comparing to the interaction Hamiltonian in the quantum mechanical picture from Eq. (1) (for a particular k value) gives us g=ω0ξ/2. Substituting for the values of ω0 and ξ from Eq. (14) gives(16)
Consider first the steady-state expectation value of the heat exchanged due to the coupling of two oscillators in contact with heat baths as given by Eq. (4) in the long time limit. This can be easily shown to reduce to(17)
The corresponding steady-state expectation value of the heat flux between the two half-spaces due to the coupling of surface modes is obtained from employing Eq. (17) in Eq. (15), making the substitutions as denoted in Eqs. (14) and (16). Taking |k|d=x, the resulting expression for the steady-state heat flux between two half-surfaces reduces to(18)
To demonstrate the dynamics of heat transfer, we take the particular case of two SiC half-spaces separated by a vacuum gap. We use the following parameters for SiC : ϵ∞=6.7, ωL=969 cm−1; ωT=793 cm−1 and Γ=4.76 cm−1 from which we obtain the surface polariton frequency for a single half-space ω0 ≈ 948 cm−1. The master equation approach can thus be used to describe dynamic evolution of the system for time scales t>>30 fs. We choose SiC as the material of our half-spaces as the surface polariton frequency falls in the infrared frequency spectrum at which it can be excited thermally at room temperature (as evidenced from the peak of blackbody spectrum at 300 K). In Figure 2A, we plot the heat flux from Eq. (15) as a function of the nondimensional time tΓ at 300 K and for vacuum gap d=5 and 100 nm. The values of heat transfer from Eq. (15) can be compared to the well-known expression for the steady-state heat transfer derived using fluctuation-dissipation theorem, PFE, which in the limit |k|d → 0 reads (19)
As seen in Figure 2A, the heat transfer rate from Eq. (15) conforms to the value predicted from Eq. (19) for time intervals tΓ 5 (see Appendix for analytical proof), and for smaller time intervals tΓ ≈ 0.1, the heat flux is observed to reach as high as ≈1.7 times the steady-state value. The heat transfer contribution as a function of in-plane wavevector component |k| at different time intervals tΓ=0.5, 2, and 6 are also shown in Figure 2B–D, respectively, and compared to the steady-state distribution obtained from using Eq. (17) in Eq. (15). In Figure 2B and C, a transient oscillatory behavior in the heat flux contribution is observed, which is indicative of the presence of a strong coupling between the surface polaritons and can be confirmed from the plot in Figure 3, where values of the ratio g(k, d)/γ>1 indicate strong coupling regime. However, the total heat flux P in Figure 2A is positive for all time scales due to the positive bias along the temperature gradient.
An approximate value of the nondimensional coupling parameter ξ=2g/ω0, relevant for heat transfer between two SiC half-spaces, can be gauged from the plot of the steady-state value of the integrand in Eq. (15) as a function of x=|k|d shown in Figure 4. The peak value at |k|d≈ 2 suggests, from Eqs. (13) and (14), that the values of ξ relevant for the coupling of two surface modes across two SiC half-spaces is approximately ≈10−2. In addition, the value of γ/ω0 ≈ 0.006 validates our assumption of weak interaction of the oscillators with the heat bath and the subsequent use of the present form of the master equation to model the dynamic effects in heat transfer .
For large wavevectors |k|>1/d, the coupling between the oscillators drops rapidly as seen in Figure 3 so that the change in thermal excitation will be negligible at such large wavevectors. This is also observed in Figure 5A and B, where a plot of 〈a†a〉/n1 and 〈b†b〉/n1 shows negligible change in the occupation numbers from the initial state at these large wavevectors when steady state is reached. In addition, in Figure 6, the occupation numbers for the two HOs are plotted as a function of time. It is observed that both oscillators attain the same steady-state occupation number within a time interval ≈10tΓ, with the temperature difference between the heat reservoirs providing the gradient for continued heat exchange in the steady state.
To conclude, we have shown that the theory of dynamics of coupled HOs connected to heat baths can be used to quantify the contribution of surface modes to the dynamic near-field heat transfer between two half-spaces separated by a vacuum gap. For two SiC half-spaces, it is observed that steady state is reached for time scales tΓ5 (which corresponds to t≈ 50 ps), and for smaller time scales, heat flux can be as high as 1.7 times the steady-state value. Experimental verification of these results would require not just the ability to measure heat transfer between macroscopic objects for small spacings but also fast response time with picosecond resolution. Recent experimental advancements show that it is possible to measure heat transfer values for gaps as low as 0.2–7 nm between an STM tip and a flat substrate  and sub-100 nm for that between flat surfaces . These measurement techniques would have to be combined with ultrafast thermometry methods (such as transient thermoreflectance technique ) to be able to measure the dynamic values of near-field heat transfer between macroscopic objects. As, with advancement in nanotechnology near-field heat transfer between objects has increasing significance, the heat transfer between components of a magnetic storage device, which are spaced a few nanometers apart is a case in point , we hope that our results, along with the other recent articles , , will pave the way for experimental verification of dynamic near-field heat transfer between objects. Although the description here has been provided for planar surfaces, the possibility of extension of this theory to other macroscopic surfaces such as microspheres and STM tips used in near-field heat transfer measurements , , ,  can also be explored.
This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement no. 702525, Funder Id: http://dx.doi.org/10.13039/501100007601. We would like to acknowledge the useful discussion with Dr. Svend-Age Biehs. The collaboration was strengthened by a 2-week visit of K.S. to the Texas A&M University. G.S.A. thanks the Welch Foundation (award no. A-1943-20180324) for supporting this work.
Here, we show the equivalence of the expression of steady-state heat transfer arrived using the coupled HO method as given in Eq. (18) to that derived from fluctuational electrodynamics principles in Eq. (19). This has been shown numerically in Figure 2 for the Lorentzian form of the dielectric function given in Eq. (10). However, the complicated expressions for the natural frequency ω′(|k|d) and the coupling constant ξ′(|k|d), as given in Eqs. (12) and (13) respectively, preclude us from showing the analytical equivalence of these two expressions for the general form of the dielectric function. For simplicity, we consider the particular case of ε∞=1 and ωT=0, where the expressions for ω′(|k|d) and ξ′(|k|d) reduce to the simple forms: and ξ′(|k|d)=e−|k|d. Such a form of dielectric function is valid for some materials such as gold (for which the corresponding parameters are ε∞=1, ωT=0, ωL=1.71×1016 s−1, Γ=1.22×1014 s−1) . We also assume that the temperature T of the half-space to be such that kBT>>ħω0. In these limits, Eq. (19) reduces to(20)
where and As the rational function M(ω)/N(ω) (which is an even function in ω) has no poles on the real axis, the integral over ω can be carried out in the complex plane using Cauchy’s residue theorem. Eq. (20) then reduces to(21)
By making the substitutions ξ′(|k|d)=e−|k|d and x=2|k|d in Eq. (18), the expression for PFE(d) in Eq. (21) can be shown to match that for Pst(d) in Eq. (18). A similar derivation for showing this equivalence can be found in Ref. .
Polder D, Van Hove M. Theory of radiative heat transfer between closely spaced bodies. Phys Rev B 1971;4:3303–14.Google Scholar
Loomis JJ, Maris HJ. Theory of heat transfer by evanescent electromagnetic waves. Phys Rev B 1994;50:18517.Google Scholar
Pendry J. Radiative exchange of heat between nanostructures. J Phys Condensed Matter 1999;11:6621.Google Scholar
Joulain K, Mulet J-P, Marquier F, Carminati R, Greffet J-J. Surface electromagnetic waves thermally excited: radiative heat transfer, coherence properties and Casimir forces revisited in the near field. Surf Sci Rep 2005;57:59–112.Google Scholar
Sasihithlu K, Narayanaswamy A. Proximity effects in radiative heat transfer. Phys Rev B 2011;83:161406.Google Scholar
Chance R, Prock A, Silbey R. Molecular fluorescence and energy transfer near interfaces. Adv Chem Phys 1978;37:65.Google Scholar
Le Ru E, Etchegoin P. . New York, USA, Elsevier, 2008.Google Scholar
Jones AC, Raschke MB. Thermal infrared near-field spectroscopy. Nano Lett 2012;12:1475–81.Google Scholar
Babuty A, Joulain K, Chapuis P-O, Greffet J-J, De Wilde Y. Blackbody spectrum revisited in the near field. Phys Rev Lett 2013;110:146103.Google Scholar
Pendry JB. Negative refraction makes a perfect lens. Phys Rev Lett 2000;85:3966.Google Scholar
Otey CR, Lau WT, Fan S. Thermal rectification through vacuum. Phys Rev Lett 2010;104:154301.Google Scholar
Van Kampen N, Nijboer B, Schram K. On the macroscopic theory of van der Waals forces. Phys Lett A 1968;26:307–8.Google Scholar
Sernelius BE. Surface modes in physics. Weinheim, Germany, Wiley-Vch, 2011.Google Scholar
Biehs S-A, Agarwal GS. Dynamical quantum theory of heat transfer between plasmonic nanosystems. J Opt Soc Am B 2013;30:700–7.Google Scholar
Barton G. Classical van der Waals heat flow between oscillators and between half-spaces. J Phys Condensed Matter 2015;27:214005.Google Scholar
Yu R, Manjavacas A, de Abajo FJG. Ultrafast radiative heat transfer. Nat Commun 2017;8:2.Google Scholar
Principi A, Lundeberg MB, Hesp NC, Tielrooij K-J, Koppens FH, Polini M. Super-Planckian electron cooling in a van der Waals stack. Phys Rev Lett 2017;118:126804.Google Scholar
Dorofeyev I. Coupled quantum oscillators within independent quantum reservoirs. Can J Phys 2013;91:537–41.Google Scholar
Ghesquière A, Dorlas T. Entanglement of a two-particle Gaussian state interacting with a heat bath. Phys Lett A 2013;377:2831–9.Google Scholar
Maier SA. Plasmonics: fundamentals and applications. Berlin, Germany, Springer Science & Business Media, 2007.Google Scholar
Agarwal GS. Quantum optics. Cambridge, UK, Cambridge University Press, 2013.Google Scholar
Walls DF. Higher order effects in the master equation for coupled systems. Z Phys A Hadrons Nuclei 1970;234:231–41.Google Scholar
Carmichael HJ. Statistical methods in quantum optics 1: master equations and Fokker-Planck equations. New York, USA, Springer, 2003.Google Scholar
Palik E. Handbook of optical constants of solids. Cambridge, MA, USA, Academic Press, 1985.Google Scholar
Kloppstech K, Könne N, Biehs S-A, et al. Giant heat transfer in the crossover regime between conduction and radiation. Nat Commun 2017;8. Article number – 14475. https://www.nature.com/articles/ncomms14475.Google Scholar
Song B, Thompson D, Fiorino A, Ganjeh Y, Reddy P, Meyhofer E. Radiative heat conductances between dielectric and metallic parallel plates with nanoscale gaps. Nature Nanotechnol 2016;11:509.Google Scholar
Caffrey AP, Hopkins PE, Klopf JM, Norris PM. Thin film non-noble transition metal thermophysical properties. Microscale Thermophys Eng 2005;9:365–77.Google Scholar
Challener W, Peng C, Itagi A, et al. Heat-assisted magnetic recording by a near-field transducer with efficient optical energy transfer. Nat Photonics 2009;3:220.Google Scholar
Rousseau E, Siria A, Jourdan G, et al. Radiative heat transfer at the nanoscale. Nat Photonics 2009;3:514–7.Google Scholar
Shen S, Narayanaswamy A, Chen G. Surface phonon polaritons mediated energy transfer between nanoscale gaps. Nano Lett 2009;9:2909–13.Google Scholar
Kim K, Song B, Fernández-Hurtado V, et al. Radiative heat transfer in the extreme near field. Nature 2015;528:387.Google Scholar
Chapuis P-O, Volz S, Henkel C, Joulain K, Greffet J-J. Effects of spatial dispersion in near-field radiative heat transfer between two parallel metallic surfaces. Phys Rev B 2008;77:035431.Google Scholar