Near-field thermal upconversion and energy transfer through a Kerr medium : Theory

We present an approach for achieving large Kerr $\chi^{(3)}$--mediated thermal energy transfer at the nanoscale that exploits a general coupled-mode description of triply resonant, four-wave mixing processes. We analyze the efficiency of thermal upconversion and energy transfer from mid- to near-infrared wavelengths in planar geometries involving two slabs supporting far-apart surface plasmon polaritons and separated by a nonlinear $\chi^{(3)}$ medium that is irradiated by externally incident light. We study multiple geometric and material configurations and different classes of interveening mediums---either bulk or nanostructured lattices of nanoparticles embedded in nonlinear materials---designed to resonantly enhance the interaction of the incident light with thermal slab resonances. We find that even when the entire system is in thermodynamic equilibrium (at room temperature) and under typical drive intensities $\sim\mathrm{W}/\mu\mathrm{m}^2$, the resulting upconversion rates can approach and even exceed thermal flux rates achieved in typical symmetric and non-equilibrium configurations of vacuum-separated slabs. The proposed nonlinear scheme could potentially be exploited to achieve thermal cooling and refrigeration at the nanoscale, and to actively control heat transfer between materials with dramatically different resonant responses.


Introduction
The field of nonlinear optics has experienced unprecedented growth in the last several decades, leading to advances in a wide range of optical technologies with applications for signal processing [1,2], detectors [3,4], spectroscopy [5], among others. Due to the inherently weak nature of bulk optical nonlinearities [6], most nonlinear devices rely on resonant systems e.g. large-etalon mirrors [7], photonic-crystal defects [8] and plasmonic resonators [9], that confine light into small mode volumes and over long timescales [10], thereby reducing power requirements. As these power requirements are scaled down, even relatively small effects stemming from thermal fluctuations can be altered by material nonlinearities [11,12], but such phenomena have only just begun to be explored [13,14]. For instance, we recently showed in Ref. [13] that at high temperatures, optical nonlinearities can alter the emission spectrum of photonic resonators, leading to assymetric lineshapes and greater-than-blackbody emission under passive (purely thermodynamic) operating conditions. In subsequent work [15], we showed that an optically driven photonic resonator can exhibit nonlinear thermal effects at a lower "temperature scale" as well as lead to new phenomena, including the appearance of Stokes/anti-Stokes emission lines and thermally activated transitions [11]. In this paper, we extend our earlier work to show that the Kerr χ (3) nonlinear response of a passive medium can be exploited to efficiently upconvert thermal radiation from mid-infrared to near-infrared or visible wavelengths. Thermal radiation from hot bodies has been studied for over a century now [16][17][18]. The recent development and application of powerful numerical techniques capable of describing thermal radiation from complex, structured materials has paved the way for the design of emitters with unique and spectrally selective properties [19][20][21]. A signficant body of work has focused on the study of radiation in the near field (short gap sizes ≪ thermal wavelength λ T ≈ 10µm near room temperature), where bound surface resonances can contribute heat and cause two objects held at different temperatures and separated by sub-micron gaps to exchange heat at rates that are orders of magnitude larger than those predicted by the Stefan-Boltzmann law (applicable only in the far field) [22][23][24][25]. This in turn has created opportunities for potential advances in the areas of photovoltaics  [33,36,37]. In this paper, we propose a different active mechanism for tailoring radiative heat transfer that exploits parametric optical nonlinearities mediated by externally incident light to extract and upconvert "thermal energy" trapped in the near field of a planar body unto another. As a proof of concept, we analyze heat exchange in planar geometries supporting surface plasmon polaritonic (SPP) resonances at far-apart wavelengths. In particular, we study planar configurations of silicon carbide (SiC) slabs separated from either silver (Ag), potassium (K), or indium tin oxide (ITO) slabs by a nonlinear (bulk or nanostructured) chalcogenide (ChG) film. Although the large SPP mismatch of the slabs leads to negligible heat exchange in the absence of a pump-SiC supports mid-infrared SPPs while the absorber slabs can support either near-infrared or visible SPPs-we show that under external illumination, energy transfer across the gap can be significant. In particular, we consider multiple interveening-gap designs-either bulk nonlinear thin films or nanostructured media consisting of lattices of nanoparticles embedded in a nonlinear medium-and find that even when the entire system is in thermodynamic equilibrium (at room temperature), the energy flux rate of mid-infrared photons which get upconverted and subsequently absorbed at near-infrared or visible wavelengths can be as large as 10 4 W/m 2 for relatively low pump intensities ∼ W/µm 2 , approaching and even exceeding typical flux rates observed in more commonly studied, nonequilibrium (passive) scenarios in which the slabs are held at a large temperature differential. Hence, this scheme allows significant thermal extraction and absorption of thermal radiation at short wavelengths (otherwise inaccessible under purely passive scenarios) and between very different (such as non-resonant) materials. Our calculations suggest that with additional design optimizations, it may be possible to exploit this approach for radiative cooling and refrigeration.
The paper is organized as follows. In Sec. 2, we consider a representative and generic system consisting of three resonators supporting modes at far-apart frequencies ω 1 , ω 2 and ω 3 = ω 1 + 2ω 2 . The resonator modes are coupled by a Kerr χ (3) medium which is driven at ω 2 by externally incident light, leading to upconversion of thermal radiation from ω 1 to ω 3 . Our analysis is based on a coupled mode theory framework [43] that describes the underlying resonant four-wave mixing process and lays out general geometric and operating conditions required to maximize upconversion. In Sec. 3, we consider concrete, physical examples of this scheme based on (a) Schematic of three equal-temperature T resonators supporting modes at ω j , each with decay rate γ j stemming from internal dissipation γ j d and/or radiation into an external channel γ jc , with j = {1, 2, 3}. The thermal modes 1 and 3 are strongly coupled to one another through a four-wave mixing process involving a Kerr χ (3) medium excited by externally incident light of power P from a waveguide that resonantly couples to the mode at ω 2 = 1 /2(ω 3 −ω 1 ). Nonlinear mixing between modes 1 and 3 can be described by an effective, linear but power-dependent coupling κ [Eq. (7)] obtained via a spatial overlap between the three modes, given in Eq. (4). (b) Thermal emission spectra P j (ω) (blue and red curves) normalized by k B T, for a choice of far-apart frequencies ω 1 = k B T/ and ω 3 = 10ω 1 , and decay rates γ j d = γ jc = 0.01ω 1 , as a function of the dimensionless frequency ω/k B T. Emission at any ω is bounded above by the Planck distribution Θ(ω, T) (black curve) in the absence of the pump κ = 0 (blue curves), but is exponentially enhanced under finite κ > 0 (red curves). (c) Thermal emission P 1 (ω) (solid lines) and heat-transfer P ex (ω) (dotted lines) spectra near ω 1 as a function of the dimensionless frequency (ω − ω 1 )/γ 1 , demonstrating splitting of the resonances into Stokes (red-shifted) and anti-Stokes (blue-shifted) peaks, which grow apart with increasing κ and further enhance emission.

Coupled-mode theory
We first illustrate the basic thermal upconversion mechanism by considering a representative system, depicted in the top inset of Fig. 1(a), involving resonators that support modes at ω 1 , ω 2 and ω 3 = ω 1 + 2ω 2 , the second of which is coupled to an external channel. The modes are assumed to have widely different frequency and therefore do not couple linearly to one another. They can however interact nonlinearly through a four-wave mixing process [6] mediated by a χ (3) medium via their field E j ( j = 1, 2, 3) overlaps and initiated by externally incident light at ω 2 from the channel. Such a system is well described by the following temporal coupled-mode equations in terms of a few key geometric parameters [43]: where a j denotes the mode amplitude of mode j ∈ [1,3], normalized so that |a j | 2 is the mode energy, and γ j = γ j d + γ jc denotes its decay rate, resulting from either material dissipation γ j d or coupling to external channels γ jc (e.g. radiation or a waveguide). Each mode is assumed to be in local thermodynamic equilibrium [13] and hence subject to thermal sources ξ j satisfying is the Planck distribution associated with the local resonator temperature T j and · · · denotes a thermodynamic, ensemble average. [39]. A monochromatic coherent drive incident from the channel, s in = s 0 exp(iω 2 t) normalized such that P = |s 0 | 2 denotes the power, facilitates the nonlinear interaction captured by the following nonlinear coupling coefficient: that depends on a spatial overlap of the linear cavity fields E αi over the (generally anisotropic) susceptibility χ (3) ijkl of the nonlinear medium in space. The indices α ∈ [1, 2, 3] and {i, j, k, l} ∈ [x, y, z] run over the mode number and cartesian components of the mode profiles, while ε 0 and ε denote the vacuum and relative permitivity of the system.
Since typical pump energies tend to be much greater than the available thermal energy in the system, i.e. P ≫ γk B T , it is safe to ignore the down-conversion term −iβω 2 a 3 a * 1 a * 2 in Eq. (2), in which case the equation for the pump decouples and is linear in the incident drive field: Such an undepleted-pump approximation [45] greatly simplifies the description of the four-wave mixing process, which is thence described by the following coupled linear equations: Here, the role of the mediator mode at ω 2 is captured by the effectively linear coupling coefficient, Generally, in addition to frequency mixing, the Kerr nonlinearity also leads to cross-phase modulation [6], which acts to shift the resonator frequencies and hence disturbs frequency matching, ω 3 = ω 1 +2ω 2 [45,46]. Since the corresponding modulation terms are temporally and hence dynamically decoupled from the effective equations under the undepleted approximation, in the following we account for such as well as other possible sources of frequency mismatch, e.g. material dispersion or even fabrication imperfections, by introducing a time-independent frequency offset ∆ω = ω 3 − ω 1 − 2ω 2 into one of the cavity frequencies.
To analyze energy transfer in this system, it suffices to consider the rate of energy loss associated with each mode, d |a j | 2 dt , obtained through Eq. (5) and Eq. (6). Collecting terms proportional to the coupling coefficient κ, one finds that the rate of energy extraction from a 1 is given by P ex = 2 Im[κ * exp(2iω 2 t)a * 3 a 1 ] , while the rate at which energy is upconverted and gained by a 3 is given by P up = 2 Im[ ω 3 ω 1 κ * exp(2iω 2 t)a * 3 a 1 ] = ω 3 ω 1 P ex . Note that in contrast to the case of two resonantly and linearly coupled modes [32], energy exchange within this four-wave mixing process is not symmetric, i.e. P ex P up , but instead satisfies a photon-number conservation condition [44], i.e. P ex ω 1 = P up ω 3 , which ensures that the number of photons lost by a 1 is equal to that gained by a 3 . Finally, the linearity of the coupled-mode equations allows the extraction/upconversion rates and emission rates P j = 2γ jc |a j | 2 to be expressed in closed form, leading to the following power spectral densities: , and as noted above, the upconverted power P up = ω 3 ω 1 P ex . The expressions above capture the most important features of this four-wave mixing scheme.
As an example, we consider the particular scenario of equal-temperature T cavities with ω 1 = k B T / , ω 3 = 10ω 1 , and zero frequency mismatch ∆ω = 0. We focus on the special situation of resonances having equal dissipative and radiation rates, γ j d = γ jc = 0.01ω 1 , in which case both resonators exhibit perfect thermal emissivities ϕ j ≡ P j /Θ(ω j , T ) (blue curves) in the absence of the pump, i.e. κ = 0. Note that we have chosen ω 3 ≫ k B T / = ω 1 , in which case there is negligible radiation at ω 3 despite the near-unity emissivity. Figure 1(b) shows the thermal radiation spectrum near the two resonances for both κ = 0 (blue curve) and κ = γ 1 (red curves), revealing giant enhancements in ϕ 3 ≫ 1 in the presence of the pump due to significant energy transfer from ω 1 to frequencies ω 3 ≫ k B T / at which fluctuations are otherwise exponentially suppressed by the Planck distribution (black curve). Figure 1(c) shows P 1 (ω) and P ex (ω) near ω 1 , illustrating that the pump also causes the mode resonances to split into Stokes (−) and anti-Stokes (+) peaks which grow apart with increasing κ and have center frequencies: As a consequence, in addition to mediating energy transfer, the pump-induced red shift allows the resonator to effectively draw additional energy available at longer wavelengths from the Planckian reservoir, thus increasing its emission rate. Specifically, owing to the red shift, the largest emissivity associated with the Stokes mode is found to be ϕ max . Such enhancements, however, are achieved only in the unrealistic limit of strong coupling κ ≫ γ 1,3 , perfect frequency-and rate-matching, ∆ω = 0 and γ 1 = γ 3 , respectively, and negligible spurious losses, γ 1c = γ 3d = 0. The largest possible extraction rate in this system, P max ex = 2γ 1d k B T , occurs under similar conditions, except in a regime wherein thermal excitations at ω 1 are upconverted and reabsorbed at ω 3 at a faster rate than they decay in resonator 1, achieved in the limit of |κ| ≫ γ 3 ≫ γ 1 . While such a strong-coupling limit typically renders the coupled-mode framework invalid [43], as we show below, it is still possible to observe significant splitting in situations κ γ j where coupled-mode theory is still accurate. Moreover, while perfect frequency matching ∆ω = 0 is desirable in order to guarantee optimal flux rates, in practice one can still achieve significant transfer rates so long as the frequency mismatch is smaller than either the resonator bandwidths, ∆ω γ j , or the coupling rate, κ γ j . Finally, as expected, it follows from Eq. (7) that in order to achieve strong coupling |κ|, one requires large nonlinear overlaps β, input power P, and long lifetimes 1/γ 2 .
The coupled-mode equations above describe a wide variety of triply resonant systems, with (a) Schematic of two plates consisting of different materials held at temperatures T 1 and T 3 and which support surface plasmon polaritons (SPPs) at far-apart resonant frequencies ω 1 and ω 3 . Also shown is a schematic of the corresponding SPP dispersions ω j (k) (blue and red). Thermal upconversion and transfer of energy from ω 1 to ω 3 is facilitated by a mediator mode at ω 2 ∼ (ω 3 − ω 1 ) 2 (green). (b) Table summarizing three different geometric configurations that result in significant thermal upconversion, along with various possible choices of emitting, absorbing, and nonlinear materials. The main difference between configurations is the choice of interveening medium, which consist of either a (i) lattice of nanoparticles embedded in the nonlinear medium or (ii) bulk nonlinear thin film, and serves as a tunable means to enhance the incident light. the choice of implementation affecting only the corresponding coupled-mode parameters. Aside from allowing different temperature reservoirs, one reason to consider a system of three physically distinct resonators is that it allows extraction of energy from one of the resonators while mitigating additional heating introduced by the coherent drive, thus paving the way for an active near-field cooling mechanism analgous to a recently proposed scheme based on gain media [31]. In what follows, we examine a set of possible, physical implementations of this upconversion scheme based on a configuration of two planar slabs supporting a plurality of resonances and which allow significant extraction of thermal energy from one slab to another.

Near-field thermal upconversion and energy transfer between slabs
Consider the planar geometry depicted schematically in Fig. 2(a) and comprising two semiinfinite materials of relative permittivities ǫ 1 and ǫ 3 which are held at temperatures T 1 and T 3 , respectively. The slabs are separated by a gap of size d that is filled with a χ (3) nonlinear medium of permittivity ǫ 2 . They support a large number of SPPs localized around their respective interfaces and characterized by their conserved in-plane momenta, k 1 and k 3 , with resonant frequencies ω 1 (k 1 ) and ω 3 (k 3 ) satisfying, where k z j = ǫ j ω 2 /c 2 − k 2 and k = |k|. Such a configuration is known to result in large thermal energy transfer when the two slabs are identical and hence support equal-frequency resonances, ω 1 = ω 3 [38]. In what follows, we consider the atypical situation of two different materials with dissimilar dispersions and far-apart SPP resonance frequencies, where there is negligible heat exchange under passive, non-equilibrium conditions but where the presence of the nonlinear medium and a (mediator) mode at ω 2 ≈ ω 3 − ω 1 2 that is excited by externally incident light can cause a large amount of thermal energy in slab 1 (emitter) to be upconverted and absorbed in slab 3 (absorber). A schematic of the frequency dispersions of such a system is shown in Fig. 2(a). As discussed above, in order for any two modes (in principles characterized by different k 1,3 ) to couple efficiently, their frequencies must satisfy the frequency-matching condition ∆ω γ 1,3 , κ described above. While material dispersion and intrinsic fabrication imperfections generally preclude this from occuring for all modes, such a condition can however be enforced by appropiate nanostructuring and/or dispersion engineering [41,42]. In principle, while it is possible to achieve frequency matching for a plurality of thermal modes by introducing and exciting multiple resonances near ω 2 , such a situation would require a more complicated analysis and is therefore out of the scope of this work.
For the sake of generality, we consider two classes of possible geometric configurations, depicted in Fig. 2(b) and differing primarily in the choice of material and interveening medium, which involve either (i) nanostructured or (ii) bulk nonlinear materials supporting SPP resonances at ω 2 . In (i), laterally incident light couples to a periodic lattice of (composite nanospheres or nanodisks) dipolar resonances while in (ii), vertically incident light couples to a low-frequency SPP through a grating. As described below, we ensure that regardless of implementation, an emitter mode at k 1 couples to a unique absorber mode at k 3 , which greatly simplifies the calculation of flux transfer by avoiding otherwise cumbersome analysis of multiply interacting degenerate modes. Such a simplification also allows us to easily compute the nonlinear coupling coefficient corresponding to each pair of modes β(k 1 , k 3 ) and to exploit the analytical expression for P ex (ω, k 1 , k 3 ) given in Eq. (8). In either configuration, the sparsity (low filling fractions) and off-resonant nature of the nanoparticles and grating allow us to ignore their impact on the translational symmetry and dispersion relations of the slab modes.
The table in Fig. 2(b) summarizes a number of potential configurations and material choices. For the sake of comparison and consistency, we choose the emitter to be SiC, whose permittivity is obtained from [53], and the interveening nonlinear medium to be ChG, whose permittivity ǫ 2 = 6.25 and χ (3) ∼ 1 × 10 −17 m 2 /V 2 are taken from various references [61-64]. For computational and conceptual convenience, we assume an isotropic χ (3) tensor with elements χ xxxx = 3 χ xxyy = 3 χ xyxy = χ (3) [6,62]. While there are many choices of possible absorber materials, e.g. K, Ag, ITO, Au, and AZO, here we focus on a select few, depending on the choice of implementation; their permittivities are taken from several references [52, 55, 57]. Finally, we fix the gap size to be d = 60nm and consider only the situation in which the entire system is at thermodynamic equilibrium, with T 1 = T 3 = 300K, for which there is zero heat exchange in the absence of the pump.

Nanoparticles lattice
In this section, we examine situations in which the interveening medium consists of a twodimensional lattice of nanoparticles embedded in ChG. The purpose of the lattice is threefold: First, it acts as a grating which allows externally incident light to excite a given (mediator) Bloch mode at ω 2 , described by the coupling coefficient of Eq. (7), which couples SPPs at ω 1 and ω 3 . Second, the particle shapes, sizes, and materials can be exploited to engineer the resonant frequency ω 2 and associated decay rates, and thereby the pump field and coupling coefficient κ, as needed. Third, the lattice provides many degrees of freedom with which to enforce "quasiphase matching" [40] over a broad range of k, allowing the β coefficient corresponding to multiple mode pairs (k 1 , k 3 ) to be optimized.
Owing to the subwavelength size of the nanoparticles, their optical response can be treated within a quasistatic, dipolar approximation [51]. They are also placed far apart from one another, mitigating many-body scattering effects and allowing the field induced by the incident drive to be conveniently expressed as a linear superposition p E 2 (x − x p , z) of the isolated particle resonances, where x and x p denote the transverse co-ordinates and center of each particle while E 2 (x − x p , z) denotes its mode profile. The p-polarized field profiles of the planar resonances are given by E jl (z)e ik j .x , with j = 1, 3 and l ∈ {x, y, z} (nonzero components alongk andẑ directions), in which case the nonlinear coupling coefficient β for a given mode pair is given by: Here, the sum is taken with respect to the index p of each particle, {i, j, k, ℓ} denote cartesian field components, and χ (3) ijkℓ is the Kerr tensor of the background medium. The SPP dispersions ω j (k) and corresponding (purely disspative) decay rates γ j d (k) are obtained by the complex-frequency solutions of Eq. (11). Because of the relatively low in-plane and volume filling fractions as well as the off-resonant response of the nanoparticle lattice, they are assumed to have a negligible effect on the planar dispersion and mode profiles.
In order to take advantage of the large density of states available in the near field, we choose geometric and material parameters that ensure large β for as wide a range of modes as possible. We restrict our analysis to modes having equal momentum k = k 1 = k 3 , in which case the integral in Eq. (12) no longer exhibits a phase factor that depends on k 1,3 . Such a simplifying assumption is further justified by the fact that although technically the nonlinearity can couple modes propagating in different directions (k 1 k 3 ), a situation that arises in the configuration of Sec. 3.2, these typically suffer from diminished mode overlaps and hence reduce the overall upconversion efficiency. Our analysis is further simplified by ensuring that the induced particle fields exhibit cylindrical symmetry about theẑ axis, in which case the field E 2 (x −x p ) and hence the coupling between each mode pair is independent of thek direction. The chosen implementions involve either a lattice of spherical nanoparticles or graphene nanodisks, respectively, which act likeẑ-oriented dipoles under illumination by light incident along either the parallel orẑ directions, respectively. These considerations allow us to write P ex (ω, k 1 , k 3 ) as P ex (ω, k), and thus to express the net heat extraction rate per unit area H across the gap as: where we further define the quantity P ex (k) = ∫ dω 2π P ex (ω, k) as the frequency-integrated flux at each k. Similarly, one can define the associated flux spectral density, Φ ex (ω) = ∫ ∞ 0 dk 2π kP ex (ω, k), by integrating instead over all wavevectors. Finally, in addition to the spatial profile of each mode pair, it is equally important to enforce the frequency-matching condition ∆ω(k) = ω 3 (k)−ω 1 (k)− 2ω 2 γ j , κ, which is generally not guaranteed and as illustrated below, depends sensitively on the choice of materials. Figure 3(a) shows a rectangular lattice of nanospheres of equal radii R and unit-cell size Λ x × Λ y that is illuminated by a z-polarized incident wave ∝ẑe −i(k 2 y+ω 2 t) , with k 2 = √ ε 2 ω 2 /c. The incident field excites dipolar modes with field profiles of the form E 2j (x − x p , z)e ik 2 y . Since for any given mode pair the nonlinear overlap integral in Eq. (12) involves a sum over all particle positions, vanishing unless the individual unit-cell contributions add constructively, we make an appropiate choice of period Λ y = π/k 2 . Furthermore, although the lack of a "Bloch phase" in thex direction means that at least in principle Λ x can be much smaller than Λ y , here we choose Fig. 3. (a) Schematic of a planar system of SiC and K slabs at thermal equilibrium (room temperature) separated by a gap of size d = 60nm that is filled with a rectangular lattice of nanospheres with unit-cell size Λ x × Λ y embedded in a χ (3) nonlinear medium (ChG). The SiO 2 /Au core/shell nanoparticles have core/shell radii of 15nm and 20nm, respectively, and the dimensions of the unit cell are Λ x = 5R and Λ y = πc/ω 2 . The SiC and K slabs support SPPs at frequencies ω 1 and ω 3 , respectively, which couple nonlinearly to dipolar particle resonances at ω 2 ∼ 1 2 (ω 3 − ω 1 ) excited by a laterally incident, monochromatic,ẑ-polarized planewave of frequency ω 2 and intensity I. (b) x y and yz cross sections of the particle resonances |E 2 | 2 , with yellow/black denoting maximum/zero amplitude. (c) Normalized slab mode-profiles at a representative wavenumber kd = 2, with the shaded region indicating the interveening medium. (d) Variations in the nonlinear coupling κ, frequency mismatch ∆ω, and frequency ω 1 with respect to kd, normalized by the corresponding disspation rate γ 3 of the K slab. (e) Frequency-integrated heat-extraction spectrum P ex (k), normalized by P 0 = 2γ 1 Θ(ω 1 , T), as a function of kd and for multiple incident intensities I. Also shown are the (f) associated spectral densitiy Φ ex (ω) and (g) net extracted and upconverted flux rates, H ex and H up , as a function of I. For comparison, (e-g) also show the heat-transfer rates associated with two vacuum-separated SiC slabs held at either 300K (light blue) or 1K (dark blue) temperature differences. a large enough Λ x = 5R in order to ignore many-body effects on the induced field. Taking advantage of the equal contribution of each unit cell to Eq. (12) and restricting our analysis to momentum-matched mode pairs k 1 = k 3 , the coupling κ(k) is given by:

Nanospheres
where ω 2 , σ abs , and γ 2d are the resonance frequency, absorption cross-section, and dissipation rate of the nanoparticle resonances, and I denotes the incident-field intensity. shells having inner and outer radii of 15nm and 20nm, respectively. These parameters yield a resonance frequency ω 2 = 8.74 × 10 14 rad/s, dissipative and radiative decay rates, γ 2d = 4.31 × 10 12 rad/s and γ 2c = 8.29 × 10 12 rad/s, respectively, and lattice parameters Λ x = 5R = 100nm and Λ y = πc/ω 2 = 431nm, resulting in volume and in-plane filling fractions of 0.01 and 0.03, respectively. The corresponding absorption cross-section σ abs = 1.2 × 10 −13 m 2 can also be obtained via a well-known dipolar analysis [50,51]. Choosing SiC as the emitter and K as the absorber means that the two interacting slab resonances have frequencies ω 1 ∼ 1.68 × 10 14 rad/s (λ 1 ∼ 11.22µm) and ω 3 ∼ 2 × 10 15 rad/s (λ 3 ∼ 0.94µm), respectively. Figure 3(b) shows the xy and yz cross-sections of the nanosphere mode profiles |E| 2 within a unit cell, where black/yellow represent zero/maximum values. The spatial mode profiles of the slab resonances are illustrated in (c) at a representative kd = 2, where the shaded region indicates the interveening medium. Figure 3(d) shows the variation in the coupling coefficient |κ| for increasing intensities I = {0.5, 1, 5, 10} W/µm 2 (from black to red), frequency mismatch ∆ω, and resonance frequency ω 1 , as a function of the dimensionless wavenumber kd. All quantities are normalized by the corresponding dissipation rate γ 3 at each k, which is the largest of the loss rates in the system. Figure 3(e) shows the ratio of P ex (k) to the thermal radiation rate of an isolated thermal resonance P 0 = 2γ 1 Θ(ω 1 , T 1 ) as a function of kd, while Fig. 3(f) shows the associated spectrum Φ ex (ω), in units of W/m 2 s, for multiple I. Finally, the integrated extraction H ex and upconversion H up rates are plotted in Fig. 3(g) as a function of I.
We now explain the most salient and important features associated with heat exchange in this system. As shown in (e), P ex (k) exhibits a small peak at kd ≈ 2 that grows and widens with increasing I, causing the overall heat transfer per unit area H ex to monotonically increase and eventually saturate as I → ∞. The associated spectrum Φ ex (ω) also shows a corresponding increase along with linewidth broadening with increasing I. These features are explained as follows: First, Eq. (8) shows that in the weak coupling regime |κ|/γ 3 ≪ 1, the flux rate H ex ∼ |κ| 2 ∼ I 2 . Second, as evident from Fig. 3(d), the bandwidth of the spectrum is primarily determined by the range of modes satisfying the frequency matching constraint, ∆ω γ j ,|κ|, with the peak occurring at the value of k which minimizes ∆ω. At small I or equivalently, κ/γ 3 ≪ 1, the range of modes satisfying ∆ω γ 3 is narrow, but higher intensities and hence increasing κ allow frequency matching to be satisfied over a wider range of k. Third, the rapid decrease in the flux rate at large kd ≫ 1 happens because β and hence κ depend on the spatial decay of the planar mode profiles, which become increasingly localized to the surface with increasing k. The precise value of k at which ∆ω(k) = 0 depends primarily on the choice of materials and geometric parameters: although its occurence at larger k would facilitate increased bandwidths due to the lack of dispersion at large k ≫ ω/c, the resulting exponential suppression in β suggests instead an optimal choice of particle parameters for a given d, here chosen at an intermediate kd ≈ 2. Note also that we only consider ∆ω arising from material dispersion and ignore power-dependent shifts due to cross-phase modulation, [45] which tend to be small and in any case can be compensated by suitable lattice parameters. Finally, as discussed above, the ratio P ex /P 0 can exceed one when the system enters the strong coupling regime γ 3 |κ| ≪ ω 1 , in which case the Stokes peak can draw additional thermal energy available at longer wavelengths. Figure 3 also shows the flux rate associated with a symmetric configuration of two identical SiC plates maintained at T 1 and T 3 and separated by the same gap size (albeit in vacuum). The net heat-transfer rate in this more typical scenario can be computed via the well-known fluctuational electrodynamics framework [23] and is given by: where r p jℓ is the p-polarized Fresnel reflection coefficient between mediums j and ℓ and where, as before, one can define and compute the frequency-integrated flux P (0) (k) and flux spectral density Φ (0) (ω). For the sake of comparison, we consider two different operating conditions corresponding to either (i) large (T 1 = 300 K, T 3 = 0 K) (light blue curves) or (ii) small (T 1 = 301 K, T 3 = 300 K) (dark blue curves) temperature differentials. From Fig. 3(e), one immediately observes that compared to the nonlinear scheme, the exponential decay in P (0) ex (k) occurs at smaller kd ≈ 1. The reasons are twofold: First, the increased proximity of the nanoparticle resonances to the slab interfaces results in slightly larger mode overlaps. Second and most importantly, increasing I and hence κ allows thermal energy in the SiC to be upconverted and then absorbed at a much faster rate than it is dissipated in the SiC. Consequently, the range of participating modes and hence the bandwidth of P ex (k) increases with increasing I. When combined with the aforementioned Stokes enhancement, the net effect is a significant increase in the nonlinear flux rates, which exceed the corresponding linear flux rates ≈ 10 4 W/m 2 at a temperature-dependent threshold intensity I c , which is approximately I c ≈ 3W/µm 2 in scenario (i) and I c ≈ 0.1W/µm 2 in scenario (ii). Finally, we note that in contrast to the typically investigated system of two vacuum-separated identical slabs, the frequency-matching wavenumber κ and hence peak value of k in the nonlinear scheme can be tuned by engineering the core/shell nanoparticle geometry or by an appropiate choice of emitter/absorber materials. Nevertheless, we find that the threshold intensities do not vary much with respect to the choice of absorber material. As shown in Fig. 3(f) (green curves), fixing the parameters of the emitter and core/shell nanoparticles but replacing K with Ag also leads to significant energy transfer (albeit slighly smaller owing to the smaller frequency-matching wavenumber).
Finally, we remark that our choice of geometry and operating conditions is by no means optimal. For instance, one could further increase the mode overlaps and frequency-matching wavenumbers by exploiting more complicated and diverse particle shapes (e.g. nanorods or even asymmetric particles) or by employing other emitter materials (e.g. hBN, Au), leading to greater efficiencies and lower power requirements. Moreover, while the choice of lattice period in this example is primarily motivated by the need to enforce quasiphase matching for a laterally incident pump, as we show in the next section, that and other geometric considerations are also highly dependent on the choice of incident direction and/or particle shape. Figure 4(a) shows a square lattice of doped graphene nanodisks of unit-cell size Λ × Λ that is illuminated by a normally incident, circularly polarized wave ∝ (x + iŷ)e −i(k 2 z+ω 2 t) , with k 2 = √ ε 2 ω 2 /c. Owing to the normal incidence, the induced fields E 2j (x −x p , z) at each lattice site have the same phase and consequently, each individual unit cell contributes equally irrespective of the lattice period. For simplicity, we choose Λ = 3R which allows us to ignore many-body scattering. As before, the cylindrical symmetry of the disks and the circular polarization of the incident light imply that β(k 1 , k 2 ) → β(k) is independent of thek direction. Our choice of graphene as opposed to other potential polaritonic/plasmonic materials (e.g. Au, Ag, etc) is motivated by the large degree of tunability in the resonance frequency ω 2 and hence frequency-matching wavenumber with respect to the Fermi energy or doping concentration of graphene, as well as by recent experiments exploring related structures [58,59]. In what follows, we choose graphene nanodisks of radius R = 20nm, fermi energy E f = 0.7eV, and intrinsic lifetime τ = 6 × 10 −13 , consistent with recent experimental measurements [59], which results in resonance frequencies ω 2 = 3.04 × 10 14 rad/s and dissipation rates γ 2d = 3.3 × 10 12 rad/s. For an incident light of intensity I, the coupling κ(k) is given by Eq. (14), where σ abs = 2 × 10 −14 m 2 is calculated from the effective dipolar susceptibility [58]. The corresponding field profiles are obtained via numerical simulation of Maxwell's equations using the boundary element method (BEM) [60], where the graphene nanodisk is assumed to have an effective thickness h = 0.5nm and dielectric

Nanodisks
ωh , with σ(ω) denoting the sheet conductivity of graphene [58]. Assuming SiC and ITO as the emitter and absorber in this configuration, respectively, one finds The graphene nanodisks have radius R = 20nm, Fermi energy E F = 0.7eV and are placed at a distance of 10nm from ITO slab. The SiC and ITO slabs support SPPs at frequencies ω 1 and ω 3 , respectively, which couple nonlinearly to nanodisk resonances at ω 2 ∼ 1 2 (ω 3 − ω 1 ) excited by a normally incident, monochromatic,x + iŷ-polarized planewave of frequency ω 2 and intensity I. (b) x y and yz cross sections of the particle resonances |E 2 | 2 , with yellow/black denoting maximum/zero amplitude. (c) Normalized slab mode-profiles at a representative wavenumber kd = 1, with the shaded region indicating the interveening medium. (d) Variations in the nonlinear coupling κ, frequency mismatch ∆ω, and frequency ω 1 with respect to kd, normalized by the corresponding disspation rate γ 3 of the ITO slab. (e) Frequency-integrated heat-extraction spectrum P ex (k), normalized by P 0 = 2γ 1 Θ(ω 1 , T), as a function of kd and for multiple incident intensities I. Also shown are the (f) associated spectral densitiy Φ ex (ω) and (g) net extracted and upconverted flux rates, H ex and H up , as a function of I. For comparison, (e-g) also show the heat-transfer rates associated with two vacuum-separated SiC slabs held at either 300K (light blue) or 1K (dark blue) temperature differences. ω 1 ∼ 1.68 × 10 14 rad/s (λ 1 ∼ 11.22µm) and ω 3 ∼ 8.4 × 10 14 rad/s (λ 3 ∼ 2.24µm). Figure 4(b) shows the xy and yz cross-sections of the nanodisk mode profiles |E| 2 within a unit cell. We remark that since the modes are largely confined to the nanodisk circumference, further close packing or potentially smaller lattice periods could potentially be employed to enhance κ. Figure 4(c) shows profiles of the planar resonances at a representative kd = 1, with the shaded region indicating the interveening medium. While most of the features observed in this configuration and illustrated in Fig. 4(d-g) are qualitatively similar to the previous implementation in Sec. 3.1.1 based on nanospheres, we emphasize some of the main differences. First, the possibility of exciting cylindrically symmetric nanoparticle resonances with normally incident light allows a reduction of the lattice period and leads to slighly reduced threshold intensities, with I c ≈ 0.5W/µm 2 in the scenario corresponding to SiC slabs held at a 300K temperature difference (light blue curves). Second, graphene offers additional (potentially dynamic) tunability with respect to the frequency-matching wavenumber, leading to large (order of magnitude) differences in the heat flux for relatively small changes in E f (not shown). Third, the two-dimensional nature of graphene nanodisks allows greater choice in their vertical placement; one could for instance consider mutliple sets of nanodisks distributed vertically throughout the interveening medium, which would further increase κ. One important constraint to consider when choosing the shapes of the nanoparticles is the choice of pump and wavelengths; for in- Fig. 5. (a) Schematic of a planar system of SiC and K slabs at thermal equilibrium (room temperature) separated by a gap of size d = 60nm that is filled with χ (3) nonlinear medium (ChG). The SiC and K slabs support SPPs at frequencies ω 1 and ω 3 , respectively, which couple nonlinearly to a SPP resonance at ω 2 ∼ 1 2 (ω 3 − ω 1 ) that is excited via a grating (of period Λ) by monochromatic light of frequency ω 2 , intensity I, and wavevector k 2y incident at an angle θ inc with respect to theẑ direction. (b) Schematic illustation of the coupling of SPP resonances of wavevectors k 1 and k 3 ; the polar plot shows the directional dependence of the nonlinear coupling coefficient κ(k, θ), where k 1 = (k, θ) is expressed in polar coordinates, with θ denoting the angle extended by k 1 with respect toŷ. (c) Normalized slab mode-profiles at a representative wavenumber kd = 2 and angle θ = 0, with the shaded region indicating the interveening medium. (d) Variations in the nonlinear coupling κ, frequency mismatch ∆ω, and frequency ω 1 with respect to kd, normalized by the corresponding disspation rate γ 3 of the K slab. (e) Frequency-integrated, angle-averaged heat-extraction spectrum P ex (k), normalized by P 0 = 2γ 1 Θ(ω 1 , T), as a function of kd and for multiple incident intensities I. Also shown are the (f) associated spectral densitiy Φ ex (ω) and (g) net extracted and upconverted flux rates, H ex and H up , as a function of I. For comparison, (e-g) show the heat-transfer rates associated with two vacuum-separated SiC slabs held at either 300K (light blue) or 1K (dark blue) temperature differences. stance, while graphene supports polaritons in the mid-infrared, other metallic nanoparticles can support resonances at higher (e.g. near-infrared or visible) wavelengths. Finally, we note that our choice of direction for the incident pump in either scenario is motivated by the desire to restrict our analysis to cylindrically symmetric nonlinear coupling coefficients (i.e. independent of the SPP propagation direction). One could also tune the lattice parameters and particle shapes so as to allow obliquely incident pumps, though this would likely lead to much more compliated β(k 1 , k 3 ). We analyze such a situation in the following section, in which we study upconversion in the case of a bulk nonlinear thin film as the interveening medium. Figure 5(a) shows a 1d y-periodic grating of period Λ resting at the interface of the ChG and K slabs. The grating is assumed to have a negligible effect on the SPPs at ω 1 and ω 3 , and chosen so that incident light of wavevector k 2y along the y direction can excite a SPP of frequency ω 2 = ω 3 (k 2y ) localized around the ChG-K interface. The angle of incidence θ inc with respect to theẑ axis is chosen so as to satisfy,

Bulk media
thus ensuring that only the first diffracted order is excited by the pump [65, 66]. The p-polarized field profiles of the planar resonances are given by E jℓ (z)e ik j .x , with j ∈ [1, 3] and ℓ ∈ {x, y, z} (nonzero components alongk andẑ), in which case the nonlinear coupling β is given by: Notably, one finds that nonzero coupling is only possible under the momentum-matching condition, k 1 + 2k 2yŷ = k 3 , represented schematically in (b). Such a condition ensures that a wave at k 1 couples to a unique wave at k 3 , in which case β(k 1 , k 3 ) → β(k, θ), where k = |k 1 | and θ is the angle extended by k 1 with respect to the y axis. Employing the translational invariance of the overlap integrals for such a process, one finds that coupling κ(k, θ) is given by: where as before, I is the intensity of the incident light and γ 2c is obtained by solving the full scattering problem [65, 66]. The net heat-transfer rate per unit area across the gap is then given by where P ex follows from Eq. (13) and the associated flux spectral density is given by In what follows, we choose a grating period Λ = 205nm in order to couple incident light at an angle θ inc = π/4 to a SPP of frequency ω 2 = 9 × 10 14 rad/s, wavenumber k 2y = 2.8ω 2 /c, and dissipation rate γ 2d = 2.3 × 10 12 rad/s. For convenience and expediency, we ignore the impact of the grating on the dispersions at ω 1,3 and assume critical coupling, γ 2c = γ 2d . Of course, any deviation from this condition would result in decreased efficiency [10]. Choosing SiC and K as the emitter and absorber materials, respectively, one finds that ω 1 ∼ 1.68 × 10 14 rad/s (λ 1 ∼ 11.22µm) and ω 3 ∼ 2 × 10 15 rad/s (λ 2 ∼ 0.94µm). Figure 5(c) shows the various mode profiles |E| 2 at a representative kd = 2 and for θ = 0, illustrating the larger spatial extent of the mediator mode owing to its effectively smaller wavenumber k 2y d ≈ 0.5. The polar plot in Fig. 5(b) shows the normalized |κ(k, θ)|/γ 3 for increasing values of kd, illustrating the directional dependence of the coupling coefficient and the lack of frequency matching stemming from material dispersion and the absence of a ω 3 (k 3 ) mode in certain directions. Figure 5(d,f,g) shows that heat exchange in this configuration is qualtiatively similar to what is observed in the case of nanoparticle lattices. Despite the similar frequency spectrum and net flux rate dependence on I, the complicated dependence of κ on both k and θ does lead to noticeably different P ex (k), shown in (e). Moreover, while this configuration has the advantage of simplicity and ease of realization, as expected the significantly decreased modal overlaps lead to much larger power requirements. For instance, the threshold intensity needed to surpass the flux rate between two SiC slabs held at a 300K temperature difference (light blue curve) is roughly an order of magnitude larger, with I ≈ 10W/µm 2 . Finally, we remark that Fig. 5(g) also show the net flux rate achieved by replacing K with Ag under the same lattice parameters and excitation conditions, indicating signficant robustness with respect to the choice of materials.

Concluding Remarks
We have presented a scheme for achieving large near-field thermal energy transfer at the nanoscale based on the nonlinear χ (3) -mediated interaction of thermal modes with externally incident light. We employed a general coupled-mode framework to predict power requirements and thermal flux rates in planar configurations consisting of different (emitter and absorber) materials separated by either a bulk or nanostructured nonlinear medium, and in which resonantly enhanced incident light at mid-infrared wavelengths upconverts mid-infrared thermal energy to either near-infrared or visible wavelengths. We find that even when the entire system is held at room temperature and at relatively low pump intensities ∼ W/µm 2 , the rate of energy upconversion can approach and even exceed 10 4 W/m 2 , which compares with typical flux rates expected in the more common situation of identical slabs separated by vacuum and held up to large ∼ 300K temperature differences. This scheme therefore not only provides a means to achieve significant energy transfer at typically inaccessible near-infrared or visible wavelengths, but also facilitates heat exchange between dissimilar materials even when these are originally in thermodynamic equilibrium. Finally, we remark that our specific choices of planar geometries and materials are by no means optimal and represent only a proof of concept. More importantly, while this scheme could potentially be exploited in nanoscale radiative cooling [31, 35, 37] and power generation [23, 29] applications, the efficiency and utility of the nonlinear process in those situations will necessarily be degrared by heating introduced by the pump and neglected in this work. Such additional heating stems from the presence of SPP resonances at the incident (midinfrared) wavelength, which necessarily cause conductive heating of the slab but can potentially be mitigated by the introduction of a small vacuum gap separating the nonlinear material from either or both of the slabs. Similarly, instead of exploiting SPP resonant enhancement of the incident light, one might also consider dielectric resonances [67, 68] in semi-transparent materials. A detailed analysis of the merits of this scheme for radiative cooling or related applications that incorporates pump-induced heating and related considerations will be the subject of future work.

Funding
This work was partially supported by the National Science Foundation under Grant no. DMR-1454836 and by the Princeton Center for Complex Materials, a MRSEC supported by NSF Grant DMR 1420541.