We modify and extend a recently proposed four-wave mixing scheme [C. Khandekar and A. Rodriguez, Opt. Express **25**(19), 23164 (2017)] for achieving near-field thermal upconversion and energy transfer, to demonstrate efficient thermal refrigeration at low intensities ∼ 10^{9}W/m^{2} over a wide range of gap sizes (from tens to hundreds of nanometers) and operational temperatures (from tens to hundreds of Kelvins). We further exploit the scheme to achieve magnitude and directional tunability of near-field heat exchange between bodies held at different temperatures.

Near field radiative heat exchange^{1–5} is important for several emerging applications and technologies, from energy conversion^{6,7} to nanoscale heat management and cooling.^{8–11} This has motivated recent efforts aimed at achieving active control of heat transfer using gain media^{12,13} or more generally, chemical potentials.^{14} Simultaneous advances in nanofabrication have also made it possible to confine light to small volumes and over long timescales,^{15,16} allowing otherwise weak optical nonlinearities to modify even low-power phenomena like thermal radiation.^{17–20} We recently proposed an alternative mechanism for controlling heat exchange^{21} that exploits nonlinear four-wave mixing to extract and upconvert “thermal energy” trapped in the near field of a planar body unto another, from mid- to near-infrared wavelengths.^{21} In particular, we showed that the combination of resonantly enhanced optical nonlinearities and large density of states associated with tightly confined surface plasmon/phonon–polariton (SPP) resonances enables high-efficiency four-wave mixing in planar materials separated by nanoscale gaps, resulting in order 10^{5}W/m^{2} upconversion rates induced by externally incident mid-infrared light of moderate intensities, on the order of 10^{12}W/m^{2}.

In this letter, we show that a similar four-wave mixing scheme can be exploited to achieve thermal refrigeration and tunable heat exchange. We begin by exploring the planar configuration shown in Fig. 1(a), comprising an emitter held at temperature *T*_{e} and supporting mid-infrared SPP resonances around frequency *ω*_{1} which is separated by a vacuum gap from an absorber held at temperature *T*_{a} and supporting near-infrared SPPs around *ω*_{3}. The absorber is coated with a thin *χ*^{(3)} nonlinear film supporting a mediator resonance at *ω*_{2} ∼ (*ω*_{3} − *ω*_{1})/2 which couples to externally incident light by way of a grating. The mediating mode facilitates resonant four-wave mixing (*ω*_{1} + 2*ω*_{2} = *ω*_{3}) between the SPP resonances, resulting in cooling of the emitter by way of upconversion and energy transfer across the gap. As shown below, in contrast to passive radiative cooling mechanisms requiring large temperature differentials *T*_{e} ≫ *T*_{a}, nonlinear upconversion allows thermal energy extraction under zero or even negative differentials (*T*_{e} < *T*_{a}), constrained only by photon-number conservation.^{21} This in turn enables thermal refrigeration, where thermal energy is made to flow from a low to high temperature body (a reversed heat engine) when the system is driven by external light, which provides the work required for energy transfer. The first part of this letter is devoted to a detailed analysis of such a refrigeration scheme, illustrating not only the various design criteria but also operating regimes needed to achieve high-efficiency refrigeration, including temperature range (*T*_{e} ∼ 10−1000 K) and gap sizes. In the second part, we extend the analysis to consider a more complicated system, depicted in Fig. 2(a), where we introduce an additional thin film on top of the nonlinear medium for the purpose of enabling appreciable heat exchange under zero external drive but finite temperature differentials *T*_{e} ≠ *T*_{a}, otherwise absent due to the large SPP frequency mismatch between the emitter and absorber. This channel can thus compete with nonlinear energy upconversion to enable tunable heat flow (in both magnitude and direction) with respect to the incident drive power.

A significant refinement in this paper with respect to our earlier work^{21} is the substitution of lossy plasmonic resonances in favor of low-loss dielectric leaky modes in the nonlinear medium. While the latter are less localized than the former, they exhibit longer radiative and absorptive lifetimes and thus result in significantly lower power requirements, on the order of 10^{9} W/m^{2} as opposed to 10^{12} W/m^{2}, while also mitigating pump-induced heating. While the choice of transparent materials around the pump wavelength *ω*_{2} mitigates heating introduced by the drive, in practice we expect that efficient thermal cooling will require a vacuum gap (not considered before^{21}) in order to further limit conductive transfer stemming from spurious heating. Our theoretical analysis is based on a coupled-mode theory framework,^{22–24} previously exploited to analyze heat transfer in linear media^{25–27} and more recently generalized to consider a broad class of weakly nonlinear resonant processes,^{28,29} that provides general operating conditions and quantitative predictions while allowing us to avoid otherwise cumbersome calculations based on nonlinear fluctuational electrodynamics.^{17} Finally, we note that our predictions extend recent work in the area of non-contact refrigeration^{12,30,31} and dynamically tunable heat exchange,^{13,32} and has analogies with more established thermoelectric cooling schemes.^{33}

## THERMAL REFRIGERATION

We first consider the planar system shown in Fig. 1(a), comprising a silica (SiO_{2}) 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_{2}S_{3}, thickness *t*, permittivity *ε*_{2} = 6.25, and isotropic Kerr coefficient *χ*^{(3)} = 10^{−17} m^{2}/V^{2}.^{37–40} (Note that we assume an isotropic Kerr coefficient, *χ*_{xxxx} = 3*χ*_{xxyy} = 3*χ*_{xyxy} = *χ*^{(3)}, 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\u2225$, 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_{2} interface of frequencies *ω*_{1a} ∼ 2 × 10^{14} rad/s and *ω*_{1b} ∼ 0.8 × 10^{14} rad/s, along with a single SPP branch localized at the AZO interface of frequency *ω*_{3} ∼ 12 × 10^{14} rad/s. Also present (not shown) is a separate mediator resonance that propagates primarily within the ChG film with frequency *ω*_{2} ≈ 5 × 10^{14} rad/s and wavevector **k**_{2} = *k*_{2}*ŷ*. This mediator mode can couple to externally incident light at *ω*_{2} and angle *θ*_{inc} by way of a thin, first-order diffraction grating of period Λ, designed to satisfy $\omega 2csin\u2061\theta inc+2\pi \Lambda =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.^{42} Furthermore, four-wave mixing between SPPs is only possible under the momentum-matching condition, **k**_{1} + 2**k**_{2} = **k**_{3},^{21} thus ensuring that only a single emitter mode at **k**_{1} couples exclusively to an absorber mode at **k**_{3}. In particular, given a set (**k**_{1}, **k**_{3}) 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* ∈ [1*a*, 1*b*, 2, 3], normalized such that |*a*_{j}|^{2} is the corresponding mode energy and *ξ*_{j} represent thermal noise sources with thermodynamic correlations $\u27e8\xi j*(\omega )\xi j(\omega \u2032)\u27e9=\Theta (\omega ,Tk)\delta (\omega \u2212\omega \u2032)$, where Θ(*ω*, *T*_{j}) = *ℏω*/[ exp(*ℏω*/*k*_{B}*T*_{j}) − 1] denotes the Planck distribution corresponding to a local bath temperature *T*_{j}.^{43} 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:^{21}

Here *I* is the drive intensity and *γ*_{2t} = *γ*_{2} + *γ*_{2c} is the overall loss rate of the mediator resonance, which includes both dissipative *γ*_{2} 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\u2212k1\u22122k2\u0177)$, allowing us to simplify *κ*_{1α}(**k**_{1}, **k**_{3}) → *κ*_{1α}(*k*, *θ*), where *k* = |**k**_{1}| and *θ* is the angle between **k**_{1} 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,^{22} leading to simple expressions for the thermal extraction $P\alpha \u21923=2\u27e8Im[\kappa 1\alpha *\u2061exp(2i\omega 2t)a3*a1\alpha ]\u27e9$ and linear heat transfer $Pa\u2192b=2\u27e8Im[\kappa l*a1aa1b*]\u27e9$ rates, along with associated power spectral densities.^{21} 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} = 300*K*. 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^{5} W/m^{2} are achievable with moderate drive intensities *I* ∼ 10^{9} W/m^{2}, illustrating over three orders of magnitude improvements in power efficiency (reduced intensity requirements) over earlier configurations^{21} 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^{33} (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 *ω*_{2} is given by $Habs=4\gamma 2\gamma 2cI/\gamma 2t2$. We choose the radiative or coupling rate of the mediator mode to be *γ*_{2c} = 10^{−5}*ω*_{2}, a very reasonable estimate based on extensive theoretical^{44–46} and experimental work on similar thin gratings.^{42} The dissipation rate *γ*_{2} ≈ Im{*ϵ*_{m}}*ω*_{2}/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^{−10} (obtained from extrapolating available data^{47}), it follows that *γ*_{2} ≪ *γ*_{2c} and as shown by red curve in Fig. 1(c), the absorbed power *H*_{abs} ≈ 4*γ*_{2}*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,^{48} 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^{−2}. 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.^{5} 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*^{2} dependence associated with linear heat exchange. Comparing *H*_{ex} in the particular case of *I* = 10^{11} W/m^{2} and equal *T*_{e} = *T*_{a} = 300 K (solid red line) to the flux rate between two SiO_{2} 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.

## TUNABLE HEAT EXCHANGE

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 *χ*^{(3)} = 10^{−18} m^{2}/V^{2} 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^{14}rad/s, obtained using the fitting procedure. In addition, the system supports a SPP branch around frequency *ω*_{3} ∼ 10 × 10^{14} 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 *ω*_{2} ≈ 4 × 10^{14} rad/s and wavevector **k**_{2} = *k*_{2}*ŷ*.

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^{12} W/m^{2}, 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^{13} W/m^{2}. 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*_{0} = 2*γ*_{1}Θ(*ω*_{1}, *T*_{e}), or the thermal power available to a single SPP, for typical values of *γ*_{1} = 4.45 × 10^{11} rad/s and *ω*_{1} = 1.78 × 10^{14} 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**_{3} = **k**_{1} + 2**k**_{2} 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*.

## CONCLUDING REMARKS

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^{9} W/m^{2}, 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^{21} 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^{14} W/m^{2} where not only coupling rates are very large but also other considerations such as optical damage threshold^{50} 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^{17} remains a challenging and interesting problem for future work.

## ACKNOWLEDGMENTS

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).