Near-field refrigeration and tunable heat exchange through four-wave mixing

We first consider the planar system shown in Fig. 1(a), comprising a silica (SiO₂) emitter separated by a vacuum gap d from an aluminum-doped zinc oxide (AZO) absorber. The associated dielectric properties are obtained from various references.^34–36 The nonlinear medium is a chalcogenide (ChG) thin film of material composition As₂S₃, thickness t, permittivity ε₂ = 6.25, and isotropic Kerr coefficient χ⁽³⁾ = 10⁻¹⁷ m²/V².^37–40 (Note that we assume an isotropic Kerr coefficient, χ_xxxx = 3χ_xxyy = 3χ_xyxy = χ⁽³⁾, purely for computational and conceptual convenience, but that more generally the nature of the relevant tensor components will depend on growth and material considerations.^38,41) The p-polarized SPP resonances in this configuration are characterized by their conserved transverse momenta k and described by mode profiles of the form $El(z)eik.x∥$⁠, where l ∈ x, y, z and x_∥ is the transverse position. Figure 1(b) shows multiple mode dispersions ω(k) arising in the above configuration, for a choice of d = 30nm, illustrating two branches of SPPs localized at the SiO₂ interface of frequencies ω_1a ∼ 2 × 10¹⁴ rad/s and ω_1b ∼ 0.8 × 10¹⁴ rad/s, along with a single SPP branch localized at the AZO interface of frequency ω₃ ∼ 12 × 10¹⁴ rad/s. Also present (not shown) is a separate mediator resonance that propagates primarily within the ChG film with frequency ω₂ ≈ 5 × 10¹⁴ rad/s and wavevector k₂ = k₂ŷ. This mediator mode can couple to externally incident light at ω₂ and angle θ_inc by way of a thin, first-order diffraction grating of period Λ, designed to satisfy $ω2csin⁡θinc+2πΛ=k2$⁠. Note that the thin size of the grating (≲ 5nm) and large frequency mismatch between the mediator and SPP resonances combine such that the grating has a negligible impact on the dispersions and resonances of the slabs.⁴² Furthermore, four-wave mixing between SPPs is only possible under the momentum-matching condition, k₁ + 2k₂ = k₃,²¹ thus ensuring that only a single emitter mode at k₁ couples exclusively to an absorber mode at k₃. In particular, given a set (k₁, k₃) of momentum-matched modes, the upconversion rates can be computed using the following coupled mode equations:

where a_j denotes the amplitude of mode j ∈ [1a, 1b, 2, 3], normalized such that |a_j|² is the corresponding mode energy and ξ_j represent thermal noise sources with thermodynamic correlations $⟨ξj*(ω)ξj(ω′)⟩=Θ(ω,Tk)δ(ω−ω′)$⁠, where Θ(ω, T_j) = ℏω/[ exp(ℏω/k_BT_j) − 1] denotes the Planck distribution corresponding to a local bath temperature T_j.⁴³ The mode frequencies ω_j(k_j) and associated decay rates γ_j are obtained from the complex, eigenfrequency solutions of Maxwell’s equations, while the nonlinear coupling coefficients κ_1α (α = a, b) describing four-wave mixing are obtained via perturbation theory^28,29 and depend on a complicated, spatial overlap of the linear profiles within the nonlinear medium:²¹ 

Here I is the drive intensity and γ_2t = γ₂ + γ_2c is the overall loss rate of the mediator resonance, which includes both dissipative γ₂ and radiative decay γ_2c (induced by the periodic grating). Note that the momentum-matching condition of nonzero coupling follows by inspection of the phase factor $ei(k3−k1−2k2ŷ)$⁠, allowing us to simplify κ_1α(k₁, k₃) → κ_1α(k, θ), where k = |k₁| and θ is the angle between k₁ and the wavevector of the mediator mode (parallel to the y axis). Also included in the coupled-mode equations is the possibility of finite linear coupling κ_l between the two emitter SPPs, which is negligible in the current configuration due to the large discrepancy between ω_1a and ω_1b but turns out to be of critical importance in the configuration of Fig. 2(b).

From these coupled-mode equations, one can find the various energy transfer rates corresponding to a given set of modes (k, θ) by considering the overall energy loss rate associated with each mode,²² leading to simple expressions for the thermal extraction $Pα→3=2⟨Im[κ1α*⁡exp(2iω2t)a3*a1α]⟩$ and linear heat transfer $Pa→b=2⟨Im[κl*a1aa1b*]⟩$ rates, along with associated power spectral densities.²¹ The net flux rates H_α→β are then given by:

where α, β ∈{a, b, 3} labels the particular flux channel.

To provide a proof-of-concept demonstration of thermal refrigeration, we first consider typical geometric and operating parameters, with d = 30nm, t = 100nm, Λ = 2μm, θ_inc = 45°, and T_a = 300K. Figure 1(c) shows the net thermal extraction rate H_ex = H_a→3 + H_b→3 corresponding to two different emitter temperatures, T_e = 300 K (blue curve) and T_e = 1000 K (black curve), as a function of drive intensity I. Evidently, large flux rates H_ex ∼ 10⁵ W/m² are achievable with moderate drive intensities I ∼ 10⁹ W/m², illustrating over three orders of magnitude improvements in power efficiency (reduced intensity requirements) over earlier configurations²¹ based on lossy plasmonic mediator resonances. Note that the transparency of the emitter at pump wavelength as well as the presence of vacuum gap and large SPP frequency mismatch imply that conductive or radiative heating of the emitter due to the pump is negligible compared to heat extraction leading to its cooling. The efficiency of such reversed heat engine (refrigeration scheme) is given by a coefficient of performance³³ (COP), defined as the ratio of thermal energy extracted to power that is lost to pump-induced heating. It follows from coupled mode equations that the absorbed pump intensity at ω₂ is given by $Habs=4γ2γ2cI/γ2t2$⁠. We choose the radiative or coupling rate of the mediator mode to be γ_2c = 10⁻⁵ω₂, a very reasonable estimate based on extensive theoretical^44–46 and experimental work on similar thin gratings.⁴² The dissipation rate γ₂ ≈ Im{ϵ_m}ω₂/2Re{ϵ_m} where ϵ_m is the complex permittivity of the nonlinear medium, is obtained from perturbation theory and agrees with the exact complex eigenfrequency solution. With Im ϵ_m ≈ 10⁻¹⁰ (obtained from extrapolating available data⁴⁷), it follows that γ₂ ≪ γ_2c and as shown by red curve in Fig. 1(c), the absorbed power H_abs ≈ 4γ₂I/γ_2c by the ultra-low loss resonance is smaller than the heat extraction rates leading to COP ≫ 1 over a varying range of intensities. While such ultra-low loss resonances have been explored extensively in different context,⁴⁸ they play an important role here in minimizing the unnecessary power dissipation and enhancing the refrigeration efficiency (COP). While ideally, COP > 1 is within reach, various non-idealities such as fabrication imperfections and spurious material loss rates may lead to effectively larger dissipation rates in actual experiments and potentially smaller values of COP ≈ 10⁻². We note that these are realistic efficiencies of all such solid-state refrigeration schemes^31,33 and are acceptable given the reliability (no moving parts) and ease of on-chip implementation in comparison to other gas based refrigeration methods.

Yet another important figure of merit is the range of operating temperatures over which it is possible to cool the emitter (independently of efficiencies). Along this vein, Fig. 1(d) shows H_ex as a function of T_e for multiple values of I, illustrating a change in the sign of the flux from positive (solid) to negative (dashed) as T_e decreases past a typical transition temperature T_e ∼ 10 K. It follows that under ideal conditions under which T_a = 300 K is held fixed, heat can be extracted from the emitter until it is cooled down to temperatures on the order of tens of Kelvin. Moreover, while we have chosen so far to focus on configurations involving very small vacuum gaps d = 30 nm, as shown in Fig. 1(e), significant flux rates can nevertheless be achieved at larger separations d ∼ 100 nm. This is an important consideration for current state of the art experiments exploring near-field heat transfer in planar geometries.⁵ Interestingly, we find that the complicated dependence of the spectral flux rate on gap and drive intensity, quantified by the coupling coefficient κ(k, θ), leads to a modified relationship between the net flux and gap size compared to the typical ∼ 1/d² dependence associated with linear heat exchange. Comparing H_ex in the particular case of I = 10¹¹ W/m² and equal T_e = T_a = 300 K (solid red line) to the flux rate between two SiO₂ plates held at a large temperature differential, T_e = 300 K and T_a = 0 K, but separated by the same gap sizes (black line), one finds not only significantly larger extraction rates but also a slower polynomial decay in the former. We finally remark that apart from the efficiency comparable to other solid-state refrigeration methods,^31,33 the design flexibility and wide range of temperature regimes achievable through this scheme could prove a viable alternative depending on the application.

We now consider a slightly modified configuration, depicted schematically in Fig. 2(a), to illustrate the possibility of exploiting four-wave mixing as a means of achieving tunable heat exchange. The modified configuration consists of a silicon carbide (SiC) emitter separated by vacuum from a aluminum-doped zinc oxide (AZO) absorber. Resting on the absorber is a composite layer consisting of an additional SiC thin film of thickness h = 20nm on top of a nonlinear gallium arsenide (GaAs) film, with nonlinear χ⁽³⁾ = 10⁻¹⁸ m²/V² and dielectric properties taken from various references.^34–36,49 The presence of the additional SiC thin film results in three branches of SPPs, out of which only two, depicted as solid lines in Fig. 2(b), correspond to SPPs localized along the SiC–vacuum interfaces. In contrast to the previous configuration, however, these SPPs have non-negligible linear coupling (κ_l ≠ 0) due to similar resonance frequencies and can therefore contribute significant linear heat exchange between the emitter and absorber. In order to incorporate both nonlinear upconversion and linear heat exchange, we obtain the unperturbed frequencies, mode profiles, and linear coupling rates κ_l of the SPPs by fitting the full linear fluctuational electrodynamics calculations to the coupled-mode equations. The dashed (solid) lines in Fig. 2(b) depict the unperturbed (perturbed) resonance frequencies, ω_1a ∼ ω_1b ∼ 1.8 × 10¹⁴rad/s, obtained using the fitting procedure. In addition, the system supports a SPP branch around frequency ω₃ ∼ 10 × 10¹⁴ rad/s that is localized along the AZO interface, shown as an inset (red line). Like before, a thin grating of period Λ = 1.78 μm is used to couple incident light at angle θ_inc = 45° to a mediator resonance in the GaAs film of frequency ω₂ ≈ 4 × 10¹⁴ rad/s and wavevector k₂ = k₂ŷ.

To provide a proof-of-concept demonstration of tunable heat exchange, we consider a typical set of parameters, corresponding to d = 50 nm, h = 20 nm, t = 100 nm, and T_e = 300 K < T_a = 400 K, which in the absence of the drive nevertheless leads to a net heat exchange across the gap directed toward the emitter. Figure 2(c) shows the net H_ex = H_a→b + H_a→3 (black line) along with individual H_a→b (blue line) and H_a→3 (red line) extraction rates, as a function of drive intensity I. Evidently, at low drive intensities I ≪ 10¹² W/m², the net extraction rate H_ex is dominated by the linear heat exchange between the SiC resonances (H_a→b < 0), becoming gradually larger due to nonlinear extraction (H_a→3 > 0) with increasing I, with the reversal in heat flow across the gap occurring at I ≳ 10¹³ W/m². Note that while intuitively one might expect a decreasing amount of linear heat flow with increasing nonlinear upconversion, we observe a non-monotonic trend in H_a→b with increasing I that demonstrates instead an increase in linear flow at low intensities. Such a non-trivial interplay between the two (linear and nonlinear) processes originates from a shift in the SiC mode frequencies that ends up enhancing the otherwise sub-optimal (due to the slight frequency mismatch) linear heat flow. While the net flux rates contain contributions from a wide set of SPPs, characterized by (k, θ), the underlying behavior of flux rates for these modes can be analyzed by inspection of the frequency-integrated spectral flux P_ex(k, θ). For illustration, Fig. 2(d) shows the angular dependence of the flux rate at a fixed kd = 0.5, with Fig. 2(e) showing the underlying angle-averaged spectrum with respect to kd at multiple drive intensities, corresponding to the points marked by circles in Fig. 2(c). For convenience, we normalize these flux rates by P₀ = 2γ₁Θ(ω₁, T_e), or the thermal power available to a single SPP, for typical values of γ₁ = 4.45 × 10¹¹ rad/s and ω₁ = 1.78 × 10¹⁴ rad/s. As illustrated in Fig. 2(d), both the magnitude and direction of the flux rate depend on the angle θ and intensity I, with the latter eventually resulting in large, positive angle-averaged flux rates at larger I. Note that there exists a range of modes, corresponding to highly acute angles (grey region), for which the momentum-matching condition k₃ = k₁ + 2k₂ can never be satisfied and for which there is no nonlinear upconversion. Finally, Fig. 2(e) shows the growing contribution and role of modes satisfying momentum-matching in the net exchange, allowing nonlinear upconversion to overwhelm the linear heat flow with increasing I.

We demonstrated a four-wave mixing scheme for active near-field heat extraction. This approach not only enables efficient nanoscale thermal refrigeration at very low temperatures ∼ 10 K and low input intensities I ∼ 10⁹ W/m², but also active control of both the magnitude and direction of heat flow across vacuum gaps. While the systems explored in this work represent only a proof of concept, we are confident that other geometries and materials could result in further improvements. We note that the coupled mode approach²¹ employed above is valid as long as the decay and coupling rates are much smaller than the resonant frequencies of SPPs. It breaks down for large intensities I ≳ 10¹⁴ W/m² where not only coupling rates are very large but also other considerations such as optical damage threshold⁵⁰ become important. While coupled mode theory circumvents the need to carry out full and repeated calculations, further analysis of such nanoscale pump-thermal mixing processes using nonlinear fluctuational electrodynamics¹⁷ remains a challenging and interesting problem for future work.

This work was partially supported by the National Science Foundation (DMR-1454836), Princeton Center for Complex Materials with funding from the NSF MRSEC program (DMR-1420541) and Cornell Center for Materials Research with funding from the NSF MRSEC program (DMR-1719875).