Nanoscale Thermal Transfer – An Invitation to Fluctuation Electrodynamics

2.1 Maxwell–Langevin Equations

Within the approximations outlined in Section 1, the radiation field in a mesoscopic medium is described by the following set of macroscopic Maxwell equations:

Here, ε=ε(r, ω) and μ describe the electric and magnetic response of the medium. The charge and current densities can be represented by (‘external’) polarisation and magnetisation fields

In the case of non-magnetic and local conductors, the key response function is the conductivity σ, and (3) and (4) apply with ε=ε₀–iσ/ω, μ=μ₀, and M=0. In some applications to thermal radiation in metals, a ‘background’ dielectric response (due to the ionic cores) is taken into account, replacing ε₀ by ε_b inside the conductor.

The boundary conditions for the fields at a smooth interface follow from the Maxwell equations (3) when understood in the sense of distribution theory: a jump in the normal component of the displacement field εE, for example, is given by the surface charge density, i.e. a δ-sheet of ρ(r), localised at the interface. Similar techniques may be applied when the response functions vary rapidly near the interface on the spatial scales relevant for the fields. One may introduce a ‘macroscopic’ model for the interface where the material equations are extrapolated from the two bulk media to the surface. The deviations from this extrapolation are called ‘excess charges’ or currents, and are represented by δ-distributions localised on the (macroscopic, idealised) interface. For details, see [15], [16].

The key idea of fluctuation electrodynamics, developed by Rytov and co-workers [1], is that the sources ρ, j are random or fluctuating fields. This is why we shall call (3) the Maxwell–Langevin equations. If the sources are non-zero on average, their values would be interpreted as ‘external’ charges as in the conventional macroscopic electrodynamics. The fluctuations around the average arise from the thermal motion of carriers in the medium and are therefore determined by its temperature and its ‘oscillator strength’ (or density of states). Within the macroscopic scheme, as explained above, it is natural to apply local thermodynamic equilibrium statistics to describe the charge and current fluctuations.

2.2 Source (Current) Fluctuations

To illustrate these concepts, let us take as an example the current fluctuations in a conductor and assume for simplicity that the mean value 〈j(r, ω)〉=0. The fluctuations of the current define its noise spectrum, as follows:

The writing 〈…〉_ω on the rhs is adopted to avoid additional notation, although it is, strictly speaking, slightly abusive. This quantity provides the (power) spectral density, i.e. the Fourier transform of the autocorrelation function (Wiener–Khintchine formula)

Note that the factor 2π in (5) is tied to this convention for the Fourier transform and that such a definition makes sense in a stationary situation only (correlations only depend on the time difference). The complex Fourier expansion used here is appropriate for ‘quantum detectors’ like a photomultiplier. Glauber’s theory of these detectors works with the Hermitian field operator E(r, t) and shows that the signal involves the convolution of E with a complex exponential e^−iωt , where ω>0 is the threshold frequency (work function). In classical electronics, the signal is convolved with a real-valued reference or circuit response. In that case, an autocorrelation function based on a symmetrised operator product is relevant. Its spectrum is symmetric and is equal to the average of its ‘quantum cousin’ in (6) over positive and negative frequencies.

Thermal equilibrium and local macroscopic electrodynamics yield the following relation between spectrum and conductivity:

where both T=T(r) and σ may vary spatially (local temperature). The real part (the conductance) of the complex admittance σ is related to ohmic dissipation, hence the name fluctuation–dissipation (FD) relation for (7). This formula is the quantum extension of Johnson–Nyquist noise; see the seminal paper by Callen and Welton on the FD theorem [34] and the review by Ginsburg [2]. The current spectrum S_j (r, ω) is plotted in Figure 2 inside a Drude conductor, σ(ω)=σ_DC/(1+iωτ). The plot covers both positive and negative frequencies: note the asymmetry of the spectrum. This is generic for operator products that are not in symmetric order and hence a typical quantum effect. Indeed, with our convention for the spectral density, quantum (or vacuum) fluctuations become visible at negative frequencies. The scheme required to detect these would be based not on absorption but on emission (spontaneous and stimulated); see, for example, [35]. The quantum current fluctuations in a Drude conductor show a broad Debye-like peak at ω~–1/τ, where τ is the electronic relaxation time. (The mean free path is v_Fτ with v_F the Fermi velocity.) Another convention is to take a symmetrised correlation function so that the spectrum becomes symmetric (grey dashed lines).

Figure 2:

Spectral density S_j (ω) of current fluctuations inside a Drude conductor (7) for positive and negative frequencies. The spectrum is normalised to its low-frequency limit at T=300 K [see (8)]; the temperature dependence of σ_DC is neglected for simplicity. The frequency is normalised to the Drude relaxation rate 1/τ (we took ℏ/τ=10 meV). The dashed lines give the spectrum for the symmetrised autocorrelation of the current.

At lower frequencies, the current noise is white and given by the Nyquist formula

It is interesting to note that for typical Drude parameters and room temperature, the ‘thermal energy–time uncertainty product’, k_BTτ/ℏ, is not far from unity. In other words, good conductors show strongly damped electrical currents at room temperature.

2.3 Field Correlations

The thermal radiation field generated by the Rytov–Langevin equations (3) depends linearly on its sources. Solutions to the homogeneous equations can be used to describe radiation incident from infinity [36]. Some authors prefer to cancel these solutions by allowing for a non-zero imaginary part in ε(r, ω) in all space (see, e.g. [37]). We follow this route and get for bodies in local equilibrium the following expression for the electric field correlation function [defined as in (5), summation over double indices]:

Here, G_kl (r, s, ω) is the electromagnetic Green tensor, which gives the electric and magnetic fields radiated by a point current source located at s:

As a function of r and the index k, G_kl satisfies the boundary conditions at macroscopic interfaces. Many different conventions for prefactors in (10) and (11) appear in the literature. Explicit expressions are available in free space [38], in a planar interface [39], [40], in multi-layer structures [41], and in spherical [42] and cylindrical [43] geometries. The book by Chew [44] provides an exhaustive discussion, including numerical techniques from computational electrodynamics. Let us simply note that in practice, it is possible to avoid the volume integration in (9) and get instead an integral over the surfaces of the bodies [9], [45], [46]. (The price to pay is the restriction to spatially constant temperatures inside the bodies.) For objects with a highly symmetric shape, an expansion into orthogonal vector field modes is possible (‘principle of sufficient symmetry’), reducing (9) to a summation over products of mode functions. It is now possible to compute Green functions in quite arbitrary geometries using finite-element techniques called ‘discrete-dipole approximation’ [47], [48], [49]. These techniques permit naturally to implement the concept of local, thermal current sources; examples of recent achievements can be found in [7], [50], [51].

From a correlation function like 〈E_k (r)E_l (r′)〉_ω , it is easy to compute the electromagnetic energy density (take k=l and r=r′), correlations for the magnetic field (take the curl with respect to both r and r′), and the Poynting vector, for example. All quantities naturally appear in the spectral domain and can be split into contributions originating from the different objects. In this way, heat currents from body A to B and back can be determined by computing the flux of the Poynting vector across a surface that separates A and B [see (20)]. The elements of the electromagnetic stress tensor [38], plotted locally, permit visualising forces (momentum transfer) between objects, providing physical insight into dispersion and Casimir forces; see [52] for examples. A seminal example of this problem is Lifshitz’ derivation of the van der Waals force between macroscopic bodies [3], [53].

2.4 Quantum Field Theory

The Maxwell–Langevin approach sketched so far provides the equations of motion of the field (operators), and relevant averages (correlation functions) are obtained from those of the source currents. The question of what would be the Hamiltonian of this theory has been raised. The quantisation scheme for the macroscopic Maxwell equations developed by the Knöll and Welsch group in Jena [37], [54] gives an answer in line with the ideas of Rytov and co-workers [1]: the information contained in the macroscopic medium response functions (dielectric function ε,...) is sufficient. We adopt here a slight rewriting of the Jena equations and keep working with the current operators j(r, ω). The basic idea is to upgrade each Fourier component to a dynamic variable. The medium+field Hamiltonian is then given by the sources alone

It is interesting that the power density of Joule absorption provides the energy (12) of medium and field. The Heisenberg equations of motion yield the seemingly trivial time evolution ~e^−iωt provided the commutator of the currents is chosen as

The relation j(r, ω)=j^*(r, –ω) is imposed as initial condition and is preserved during the evolution. The local thermal equilibrium ensemble is generated by a density operator proportional to the exponential of H_MF, taking the local temperature field T(r) under the spatial integral in (12). The FD relation with its correlation function (7) then follows.

The observable electric field is given by the source current convolved with the Green function (10), and this relation can be understood as a generalised mode expansion. In this way, the Maxwell–Langevin equations are satisfied. The coupling to external sources (atoms or molecules) proceeds in the usual way by adding interaction terms to H_MF. A multipolar coupling scheme is quite natural; for details, see the book by Buhmann [54]. An important consistency check is to recover the Pauli–Jordan commutator between the electric and magnetic fields:

This is achieved by using (13) for the source currents, the Kramers–Kronig relations for the Green function, and the following identity [55], [56]:

which follows from the classical macroscopic Maxwell equations. The frequency arguments were suppressed for simplicity, and the negative real part appearing on the rhs is due to the convention adopted here for the Green tensor [compared to the energy conservation law (19) below].

2.5 Discussion

Let us finish this basic introduction with a remark on physical interpretation. A cursory scan through the literature yields different answers as to the actual status of the current fluctuations appearing in (3) and (7). One may consider them as ‘mathematical artefacts’ whose only merit is to reproduce, in thermal equilibrium, the spectrum of the radiation field [57], similar to thermostatting Langevin forces in molecular dynamics. It has been checked that equilibrium is indeed recovered if all bodies are at the same temperature: the field then follows Planck’s blackbody spectrum, modified by the bodies’ emissivity and enhanced significantly when near-field components become detectable at subwavelength distances [57]. Right from its inception, however, fluctuation electrodynamics have been applied to non-equilibrium settings, too [4], where this reference situation is no longer available. One is then inclined to adopt the viewpoint that Rytov’s fluctuating currents are a model for ‘actual thermal fluctuations’ of matter observables that are relevant as electromagnetic sources – an approximate model, of course, as in the formulation presented above, the local macroscopic theory has been applied. However, the model has some reasonable footing given the large number of microscopic constituents that make up a condensed matter system on a mesoscopic length scale of, say, a few nanometres [see (2)]. The key assumption is, of course, that on this scale, the concept of a local temperature T(r) makes sense, and that the electromagnetically relevant quantities have fluctuations that are equilibrated at this temperature. Such a model is eventually required in any description, as described insightfully by Barton [58]:

One could, of course, try to pursue the further question of how the Langevin forces themselves might be kept functioning as envisaged; to explore this would then entail hypothesising some thermostats controlling the forces, and calculating the forces instead of making assumptions about them a priori; and so on, potentially ad infinitum. In the end, one must necessarily settle for control at some level through thermostats that are external to the system, in the sense that they impose temperatures by fiat, through dynamics that do not enter the calculation and are not spelled out.

The following formula may illustrate the problem of local equilibrium:

This sum (sometimes called the ‘Rytov split’) gives the current as the sum of Ohm’s law (‘induced current’), while the second term gives the fluctuations around this value (‘fluctuating current’, sometimes marked with subscript fl). In a non-equilibrium setting like in heat transfer, the field E(r, ω) is, in general, not in a state of local equilibrium because it contains radiation from distant sources at different temperatures. The current is then also not in local equilibrium, because its induced part is simply the linear (and local) response of the medium to this non-equilibrium field. However, the fluctuations around it (fluctuating current) are characterised by the local temperature. This has been recognised as the key assumption of fluctuation electrodynamics, for example by Barton [58], [59] who comments in [58] on Loomis and Maris [60]:

[...] who do not assume local thermal equilibrium, but only that in each half-space the noise is appropriate to the temperature of its thermostat.

It is, of course, another problem how local thermal equilibrium is established in the material, and whether this happens fast enough to talk of a single temperature for all frequencies of the local current noise. A practically minded engineer might want to use an equation of motion for the temperature distribution T(r) using a combination of heat conduction and heating/cooling by absorption/emission of thermal radiation. It seems quite obvious that such a ‘self-consistent’ model requires a separation of time scales between the ‘fast’ electromagnetic fluctuations and the ‘slow’ evolution of the temperature field. Needless to say that for many experiments on near-field radiative transfer, it is sufficient to assume a stationary temperature profile, which may even be spatially flat (across each body) if the thermal conductivity is large enough. For a counter-example on short time scales, see the pump–probe experiments quoted above [33].

As a quick estimate, let us consider a typical nanoscale heat flux that aims at competing with the solar constant, q=1360 W/m². Thermal conductivities of condensed matter are in the range κ=1…10² W/Km at room temperature. On the spatial scale L of a microdevice, one thus estimates a temperature difference

which is probably below the precision of typical thermometry. In a nanojunction, the temperature profile is thus ‘step-like’ (see Fig. 1). The radiation field in the gap between the bodies has no well-defined temperature, of course. This is similar to radiation in the upper atmosphere, which is a superposition of solar and terrestrial sources.

3.1 Typical Features on the Nanoscale

Thermal radiation is a classical topic in astrophysics that determines the luminosity and the inner structure of stars. The cosmic microwave background also provides an example where the blackbody spectrum is a very accurate description.

As one works with condensed matter sources, deviations from the Planck spectrum become apparent. Some of these are well known and are described by the classical Kirchhoff concepts of emissivity and absorption (‘grey bodies’) [61]. As one enters the near field, the changes in the spectrum become dramatic because surface resonances may appear in the range of thermal frequencies. Polar materials provide a generic example where optical phonons hybridise with the electromagnetic field into surface phonon polaritons. These resonances enhance the electromagnetic energy density, and they also polarise the radiation field. An example is shown in Figure 3 for a planar body: at the wavelength 11.36 μm, silicon carbide has surface resonances (because the permittivity is negative), which show up as a prominent ‘shoulder’ at distances d≈λ/10. At even shorter distances, the energy density follows a law 1/d², which can be explained by the electrostatic fields of fluctuating dipole moments. The surface resonance appears as an enhanced amplitude of this power tail.

Figure 3:

Electric energy density spectrum as a function of distance (top panel) and degree of polarisation (bottom panel). Parameters for SiC: ε(11.36 μm)=–7.6+0.4 i and ε(9.1 μm)=1.8+4.0 i; for glass ε(500 nm)=2.25+10⁻³ i. The energy densities follow a power law 1/d² at short distances. They have been normalised to the Planck spectrum. Figure adapted from ([11], fig. 8).

The polarisation of the near field has been quantified in the plot by the following degree of polarisation:

where ε_zz (ε_||) is proportional to the spectrum of the electric field perpendicular (parallel) to the surface, and by symmetry 〈Ex2(r)〉ω=〈Ey2(r)〉ω. Note the strong polarisation when the surface resonance is excited. In the limits of short and large distances, a partial polarisation P≈1/4 arises, but for different reasons: at short distances, it follows from electrostatics that the normal field (E_z ) has a spectrum twice as large as the parallel field. At large distances, the polarisation is due to reflection and emission at glancing angles from the surface, similar to the Brewster effect. This partial polarisation is actually an artefact of assuming a planar source of infinite extent. For a discussion of the thermal emission of a spherical source, see, for example [64], [65].

Let us mention that to measure the energy density, a local measurement with a small pick-up antenna is needed. This can be achieved with sharp tips that scatter the near field [27]. The scattered signal is dominated by the immediate environment of the local probe: a lateral resolution comparable to the distance d would be typical. In this way, the thermal radiation may provide a scanning image of local electromagnetic properties and resolve metallic nano-objects on a dielectric substrate, for example.

3.2 Energy Balance and Radiative Heat Transfer

The radiative heat flux between two objects is described in electromagnetism by the Poynting vector. In our context, it has a natural spectral representation and involves the Planck spectrum as an essential factor, due to the assumption of local thermal equilibrium of the radiation sources. The symmetry of photons propagating to the left and the right is broken by the difference in temperature. As the bodies are approached to subwavelength distances, the heat flux is enhanced due to near-field coupling.

Let us start with the energy conservation law for the Maxwell–Langevin equations. In a stationary situation (the energy density is constant), one gets

for each frequency ω. The first term on the lhs involves the radiative emission spectrum (average Poynting vector 〈S〉_ω ); the second term, the absorption in the medium. (With our sign convention, passive media have ω Im ε(ω)≥0.) On the rhs, we find the power transferred to the fields by the mechanical motion of the sources [38]. Equation (19) has been commented upon in a discussion of energy conservation in a dissipative system [66]: the electromagnetic energy lost by absorption in the medium is ‘replenished’ by the sources in the same medium, as long as local equilibrium holds. This observation illustrates the general idea of treating dissipative quantum systems with Langevin dynamics: the fluctuating forces maintain the system energy at its equilibrium value. By the same token, they prevent the commutators of the system observables from decaying.

In a typical setting of heat transfer between two bodies, we may evaluate (19) in a vacuum gap between them and get ∇·〈S〉_ω =0. In a planar geometry, this means that the normal flux is constant. We quote here the formula that results from evaluating the mixed electric and magnetic correlation function between two planar bodies, ‘hot’ and ‘cold’, at distance d [60], [67]:

where Φ(ω, T) is the flux of a blackbody. Two polarisations p are summed over, and the quantity

is proportional to the absorption of a wave incident on body b=h, c where it is reflected with amplitude r_bp=r_bp(k, ω). Finally, kz=(ω/c)2−k2 is the normal component of the wavevector in vacuum. For k>ω/c, it is purely imaginary (evanescent wave, Im k_z ≥0), the heat then transports by photon tunnelling from the bulk of one body to the other. The denominator 1−rhprcpe2ikzd in (20) describes the multiple reflection of waves between the two interfaces; its zeros correspond to ‘cavity modes’ and surface resonances. In practice, due to imperfect reflections, the modes are broadened but are still visible in the spectrum. See Figure 4 for an example, taken from the seminal paper by Polder and van Hove [67].

Figure 4:

Spectrum of heat flux between two metallic slabs, contribution of propagating waves [0≤k≤ω/c in (20)]. Upper pair of curves: electric field in the plane of incidence (p- or TM-polarisation), lower pair: electric field perpendicular to plane of incidence (s or TE). Solid lines: distance d=10 μm, dashed lines: d→∞. The peaks arise when d is an integer multiple of λ/2 (standing waves). The bodies are at temperatures T_h =315 K=T_c +1 K, and have identical conductivities described by a Drude model with σ_DC=3.9 MS/m and damping time τ=7.1 fs. Figure taken from [67], fig. 4. ©American Physical Society. Reprinted with permission.

The enhancement of the heat flux in the near field is mainly due to the contribution of evanescent waves (photon tunnelling). Taking the limit of purely non-absorbing bodies, one finds that the k-integral in (20) is limited to k≤nω/c, where n is the refractive index: in this case, the Planck spectrum (Stefan–Boltzmann law) still holds with a modified prefactor (angle-averaged emissivity). Figure 4 illustrates that the propagating modes that dominate in this range mainly redistribute the emission spectrum by forming cavity resonances. The qualitatively new power laws seen in Figure 3 emerge from deeply evanescent waves, k≫ω/c. These correspond, on the one hand, to surface polaritons with typically k~1/d where d is the distance, and on the other, to diffusion ‘modes’, k~(μ₀σ_DCω)^1/2, where (μ₀σ_DC)⁻¹ is the magnetic diffusion constant in a conductor. From an analysis of the two polarisations, these contributions have been pointed out already by Polder and van Hove [67].

4.1 Casimir Force and Conductivity

The first anomaly has been dubbed the ‘plasma vs. Drude controversy’ in the community working on dispersion forces and the Casimir effect [68], [69]. We focus here on an experiment, performed by Chen et al. [70], that does not seem to have received much attention in this discussion. In the setup, an atomic force microscope is measuring the force between a sphere (diameter ~200 μm) and a planar body. The body is a semiconductor membrane (silicon) and gets irradiated with a train of laser pulses that increases the carrier density. The density is practically stationary over the ~5-ms duration of the pulses, and the force is measured with a lock-in amplifier at the repetition rate. From the viewpoint of macroscopic electrodynamics, the irradiation changes the conductivity by a few orders of magnitude from quite small (intrinsic, related to defects, ~10 S/m) to high (~1 MS/m), while the temperature increase is negligible.

The calculated electromagnetic force agrees with the experimental data within 10%, but only for the laser-doped material. A statistically significant difference to the theory (about 50% at distance d≈100 nm) is observed for the intrinsic Si membrane. The deviation can be pinpointed to the behaviour of the material conductivity across the thermal spectrum (see [71], [72] for a discussion). Technically speaking, frequency integrals are performed as summations over the imaginary axis (Matsubara sum), using the analytical continuation of response functions. The anomaly arises from the zero-frequency term in the sum. The current discussion about this thermal Casimir anomaly has led to experiments with magnetic materials (nickel) where the deviation to theory amounts to a few orders of magnitude [73]. An explanation for this anomaly would turn precision measurements of forces in the submicrometre range into tighter constraints on gravity on small scales (including extensions of the standard model like the ‘fifth force’ and additional dimensions) [74], [75], [76].

4.2 Heat Transfer Below 10 nm

Two recent experiments on heat transfer between sharp, metallised tips and a substrate, held at different temperatures, suggest a similar anomaly when compared to fluctuation electrodynamics [77], [78]. The setups differ mainly in the size of the tips (radius of curvature). The group of Reddy and Meyhofer found good agreement with theoretical calculations done by Kim et al. [77], within uncertainties related to surface roughness. Distances down to a few nanometres were reached (tip radius ≈450 nm). Kloppstech et al. observed a heat flux much larger than theory (computed by Biehs and Rodriguez) [78], using a sharper tip (radius ≈30 nm) in a similar range of distances. The data show an onset of the ‘giant heat flux’ at d≈5…6 nm, with a roughly linear increase as d is reduced. According to calculations, a significant fraction of the heat flux should originate from the ‘shaft’ of the tip (conical shape with ~300 nm height). It has been checked that the difference between experiment and theory is not reduced when a non-local conductivity (spatial dispersion) is applied [78]. The solution to this anomaly in short-range heat transfer is currently under investigation [79], [80]. Molecular dynamics simulations like those performed by Domingues et al. [81] and Chiloyan et al. [82] are likely to provide a versatile tool to ‘bridge the gap’ between the mesoscopic and microscopic scales and to include, for example, heat transport from phonons [83].