Dynamic near-field heat transfer between macroscopic surfaces

The dynamic heat transfer between two half-spaces separated by a vacuum gap due to coupling of their surface modes is modelled using the theory that describes the dynamic energy transfer between two coupled harmonic oscillators each separately connected to a heat bath and with the heat baths maintained at different temperatures. The theory is applied for the case when the two surfaces are made up of a polar crystal which supports surface polaritons that can be excited at room temperature and the predicted heat transfer is compared with the steady state heat transfer value calculated from standard fluctuational electrodynamics theory. It is observed that for small time intervals the value of heat flux can reach as high as 1.5 times that of steady state value.

The dynamic heat transfer between two half-spaces separated by a vacuum gap due to coupling of their surface modes is modelled using the theory that describes the dynamic energy transfer between two coupled harmonic oscillators each separately connected to a heat bath and with the heat baths maintained at different temperatures. The theory is applied for the case when the two surfaces are made up of a polar crystal which supports surface polaritons that can be excited at room temperature and the predicted heat transfer is compared with the steady state heat transfer value calculated from standard fluctuational electrodynamics theory. It is observed that for small time intervals the value of heat flux can be significantly higher than that of steady state value. 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 [1][2][3][4]. This heat transfer comprises of contributions from both long-range radiative modes as well as near-field evanescent and surface modes [4,5]. 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 [4]. 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 infra-red wavelength around 10 µm and can thus be thermally excited at these temperatures.
In general, resonant excitation of surface modes plays an important role in several phenomena and applications including: decreased lifetime of molecules close to metal surfaces [6], surface enhanced raman spectroscopy [7], thermal near-field spectroscopy [8,9], concept of perfect lens [10] and thermal rectification [11]. The study of coupling of surface modes across surfaces is significant as it not only plays an important role in heat transfer but also in the van der Waals and Casimir force between them [12,13]. A coupled harmonic oscillator description for the heat transfer between nanoparticles due to coupling of surface modes was arrived at by Biehs and Agarwal [14] and estimates for both dynamic and steady state heat transfer values were arrived at. Barton [15] has considered the heat flow between two harmonic oscillators 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., [16] have recently analysed 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 resultant 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., [17] where they observe thermalization in picosecond timescales. In this paper we model the dynamic heat transfer contribution from coupling of surface modes across two dielectric planar surfaces using the master equation description of two coupled-harmonic oscillators interacting with their respective heat baths and compare the steady-state results with those obtained from fluctuational electrodynamics principles available in literature [2,4]. Our work differs from that in Ref. [16,17] in that we are interested in analysing 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. Fig. 1. Coupling between two surface modes each connected to its own heat bath that is analysed in this work. The two half-spaces are assumed to be made of the same material which is local and dispersive so that the complex dielectric function varies only with the frequency ω. Only modes (denoted by the amplitudes a(k) and b(k)) of the same in-plane wavevector component k couple across the vacuum gap d. The coupling constant g between the surface modes is derived using Maxwell's equations The paper is arranged as follows: in Section 1 the results of the dynamics of heat transfer between two coupled HO each in contact with a heat bath are summarised. While the theory has been described in detail in Ref. [14], for sake of completion the main results have been reproduced here with added details. Such a system has also been analysed previously in the context of analysing the dependence of mean interaction energy on the temperatures of heat baths [18], and entanglement between two particles [19]. In Section 2 the theory of dynamics of heat transfer between two coupled HO is extended to that between two half spaces by analysing the coupling between two interacting planar surface modes using Maxwell's equations. In Section 3 numerical values for the heat flux derived in Sec. 1 and 2 are plotted for the particular case of two silicon carbide half-spaces and compared with calculations from fluctuational electrodynamics principles.

Heat transfer dynamics between two harmonic oscillators
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 + ik z z), where the in-plane component k = (k x , k y ) and the z-component k z of the wavevector are related as k 2 + k 2 z = (ω/c) 2 with c being the velocity of light, ω the frequency of the planar wave in the vacuum gap of this system, and k being limited by |k| > ω/c. The surface polariton exists for the (k, ω) pair that satisfies the well known dispersion relation [20] |k| = (ω/c) ε(ω)/(ε(ω) + 1) so that the frequency ω 0 of the surface polariton is characterized by the wavevector k. Thus one can characterize the surface excitations in terms of oscillators with frequency ω 0 (k), complex amplitude a(k) and half line width γ(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 Fig. 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 the wavevector k is conserved for fields across planar interfaces. This implies an effective coupling of the form: where b(k) is the complex amplitude of the surface polariton on the half-space B, g(k) is the coupling constant (units of rad s −1 ) and we neglect non-resonant contributions to the coupling (i.e. we use the rotating wave approximation). Thus the mode a(k) couples to the mode b(k) only. There are no coupling terms of the form a † (k 1 )b(k 2 ), k 1 = k 2 .
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 non-retarded limit the surface modes occur for all values of k. This coupling results in a heat flow due to interaction between the modes which is a dynamical process, and which 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 [21] and is given by: Here ρ S is the density matrix for the two oscillator system. The operators a, a † , b, 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. The half-space A is at temperature T so that n 1 is the excitation of the mode a i.e., The 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 [22]. 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 [23]. From this, and relating expectation value of an operator G to the density matrix ρ S using the standard relation ∂ G /∂t = Tr(∂ρ S /∂t G) we obtain the rate of change of the excitation of surface polariton in B to be: 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 (in units W) to surface polariton in B as: To find P (t) we follow the procedure used to obtain Eq. 3 to get the dynamic matrix equation: 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 = n 1 ; 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 analyse the dynamical evolution of the system. We thus get the expression for heat transfer from Eq. 4 as:

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 line-width γ, coupling constant ξω 2 0 , and displacements x 1 (ω), x 2 (ω) where x 1/2 (ω) = e iωt x 1/2 (t) dt which are related as: The eigenfrequencies for such a system are are given by: To find the 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 with that of Eq. 8 for two coupled harmonic oscillators. Consider two half-spaces separated by a gap of width d occupying the regions z < −d/2 and z > d/2 as shown in Fig. 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| → ∞ [20]. In this electrostatic limit, satisfying continuity of perpendicular component of displacement field gives the condition for surface modes as: where ε(ω) is the dielectric function of the half-space. Taking the Lorentz model for the dielectric function: 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: where, and In the limit of large gaps |k| d → ∞ , ω 0 (k, d) ≈ ω 0 where ω 0 is the surface polariton frequency for a single half-space given by ω 0 = (ε ∞ ω 2 L + ω 2 T )/(ε ∞ + 1) and the coupling parameter ξ (k, d) ≈ 0. For the case when ε ∞ = 1 these expressions reduce to the simple forms: ω 0 = (ω 2 L + ω 2 T )/2 and ξ (k, d) = e −|k|d (ω 2 L − ω 2 T )/(ω 2 L + ω 2 T ) [15]. 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 [12] and the steady state van der Waals heat transfer [15]. By comparing Eq. 11 with Eq. 8 we can relate parameters of the harmonic oscillator model with those of the surface modes: We can thus use the results from Sec. 1 to find the heat transfer between the two flat surfaces due to coupling between the surface modes. As noted earlier, we add the contributions from all the harmonic oscillator modes (labelled 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) (Wm −2 ) is then given by: where P (t, k, d) is the rate of heat transfer (units of W) between two harmonic oscillators 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 harmonic oscillator system is given by: where the interaction term is given from: H I = ω 2 0 ξx 1 x 2 . The equivalent quantum mechanical model can be got by substituting the displacements x 1 and x 2 in terms of the commutator operators a,a † , b,b † using the standard relations x 1 = (a + a † ) /(2ω 0 ), and x 2 = (b + b † ) /(2ω 0 ). Including only terms resulting in photon exchange we get: Comparing with 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: 3. Results Consider first the steady state expectation value of the heat exchanged due to 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: The corresponding steady state expectation value of heat flux between the two half-spaces due to coupling of surface modes is obtained from employing Eq. 17 in Eq. 15, making the substitutions as denoted in Eq. 14 and Eq. 16. Taking |k|d = x the resulting expression for the steady state heat flux between two half surfaces reduces to: For demonstrating 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 [24]: ∞ = 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 Fig. 2a we plot the heat flux from Eq. 15 as a function of the non-dimensional time tΓ at 300 K and for vacuum gap d = 5 nm and 100 nm. The values of heat transfer from Eq. 15 can be compared with the well-known expression for the steady-state heat transfer derived using fluctuation-dissipation theorem, P FE , which, in the limit |k|d → 0, reads [2]: As seen in Fig. 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 Fig. 2b, 2c and 2d respectively and compared with the steady state distribution obtained from using Eq. 17 in Eq. 15. From Fig. 2b and 2c a transient oscillatory behavior in the heat flux contribution is observed which is indicative of the presence of strongcoupling between the surface polaritons, and which can be confirmed from the plot in Fig. 3 where values of the ratio g(k, d)/γ > 1 indicates the strong-coupling regime. However, the total heat flux P in Fig. 2a is positive for all time scales due to the positive bias along the temperature gradient. for 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 [22].
For large wavevectors |k| > 1/d the coupling between the oscillators drops rapidly as seen in Fig. 3 so that the change in thermal excitation will be negligible at such large wavevectors. This is also observed in Fig. 5(a) and (b) where a plot of a † a /n 1 and b † b /n 1 shows negligible change in the occupation numbers from the initial state at these large wavevectors when steady state is reached. In addition, in Fig. 6 the occupation numbers for the two  harmonic oscillators 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 harmonic oscillators connected to heat baths can be used to quantify the contribution of surface modes to the dynamic near-field heat transfer between two halfspaces 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 picoseconds) 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 a STM tip and a flat substrate [25], and sub 100 nm for that between flat surfaces [26]. These measurement techniques would have to be combined with ultrafast thermometry methods (such as transient thermoreflectance technique [27]) in order to be able to measure the dynamic values of near-field heat transfer between macroscopic objects. Since, 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 few nanometers apart is a case in point [28] -we hope that our results, along with the other recent articles [16,17], will pave the way for experimental verification of dynamic near-field heat transfer between objects. While the description here has been provided for planar surfaces, possibility of extension of this theory to other macroscopic surfaces such as microspheres and STM tips used in near-field heat transfer measurements [25,[29][30][31] can also be explored.

Appendix
Here, we show the equivalence of the expression of steady-state heat transfer arrived using the coupled harmonic oscillator method as given in Eq. 18 with that derived from fluctuational electrodynamics principles -Eq. 19. This has been shown numerically in Fig. 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 Eq. 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: ω (|k|d) = ω L / √ 2 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 × 10 16 s −1 , Γ = 1.22 × 10 14 s −1 ) [32]. We also assume that the temperature T of the half-space to be such that k B T ω 0 . In these limits Eq. 19 reduces to: where, M (ω) = Γ 2 ω 2 L ω 2 and N (ω) = 16Γ 4 ω 4 − 4ω 4 − 4ω 2 L ω 2 + (1 − e −x )ω 4 L 2 +8Γ 2 ω 2 4ω 4 − 4ω 2 L ω 2 + (1 + e −x )ω 4 L . Since 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. Equation 20 then reduces to: By making the substitutions ω (|k|d) = ω L / √ 2, ξ (|k|d) = e −|k|d and x = 2|k|d in Eq. 18 the expression for P FE (d) in Eq. 21 can be shown to match that for P st (d) in Eq. 18. A similar derivation for showing this equivalence can be found in Ref. [15].