## Abstract

Nearly all thermal radiation phenomena involving materials with linear response can be accurately described via semi-classical theories of light. Here, we go beyond these traditional paradigms to study a *nonlinear* system that, as we show, requires quantum theory of damping. Specifically, we analyze thermal radiation from a resonant system containing a *χ*^{(2)} nonlinear medium and supporting resonances at frequencies *ω*_{1} and *ω*_{2} ≈ 2*ω*_{1}, where both resonators are driven only by intrinsic thermal fluctuations. Within our quantum formalism, we reveal new possibilities for shaping the thermal radiation. We show that the resonantly enhanced nonlinear interaction allows frequency-selective enhancement of thermal emission through upconversion, surpassing the well-known blackbody limits associated with linear media. Surprisingly, we also find that the emitted thermal light exhibits non-trivial statistics (*g*^{(2)}(0) ≠ ~2) and biphoton intensity correlations (at two *distinct* frequencies). We highlight that these features can be observed in the near future by heating a properly designed nonlinear system, without the need for any external signal. Our work motivates new interdisciplinary inquiries combining the fields of nonlinear photonics, quantum optics and thermal science.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Nonlinear interactions between light fields are typically weak inside bulk media but they can be significantly enhanced in resonant systems to observe them at low optical powers. As the power requirements are scaled down, the nonlinearities can influence not only quantum optical processes [1] but also low power thermal radiation phenomena (light fields generated by thermal fluctuations of charges). The nonlinear mixing of thermal photons is a theoretically-challenging, largely-unexplored topic which is relevant in the field of thermal radiation with applications for renewable and energy-conversion technologies [2,3]. Recently, it was pointed out using semi-classical Langevin theory that resonantly enhanced Kerr $\chi ^{(3)}$ nonlinearities can be harnessed for spectral-engineering thermal radiation [4,5] and near-field radiative heat transport [6,7]. In that context, the problem of thermal radiation from a *passive* system of coupled resonators of *distinct* frequencies remained unsolved. We solve it here using quantum theory for a system containing a $\chi ^{(2)}$ nonlinear medium and discover new fundamental possibilities for the field of thermal radiation.

Our system depicted in Fig. 1 consists of two resonators at frequencies $\omega _1$ and $\omega _2 \approx 2\omega _1$, and contains a $\chi ^{(2)}$ nonlinear material. In stark contrast to prior works in the quantum-optics literature where an external drive is used [8–10], here the resonators are driven by low power thermal fluctuations inside the medium. Because of $\chi ^{(2)}$ nonlinearity, two photons at $\omega _1$ can get upconverted to yield a single photon at $\omega _2$, and a single photon at $\omega _2$ can get downconverted to yield two photons at $\omega _1$. The resonantly enhanced, thermally driven upconversion and downconversion processes can alter the well-known Planck’s distribution as illustrated in Fig. 1, causing redistribution of thermal energy in different parts of the frequency spectrum.

It turns out that the analysis of thermal radiation from such nonlinearly coupled resonators of *distinct* frequencies is not trivial. Theoretically, it requires careful examination of their thermal equilibrium behavior. It is known that the individual, uncoupled resonators (harmonic oscillators) in equilibrium with a reservoir at temperature $T$ contain an average thermal energy given by the Planck’s function $\Theta (\omega ,T)=\hbar \omega /[\textrm {exp}(\hbar \omega /k_B T)-1]$, above the vacuum zero point energy [11]. In the classical regime ($\hbar \omega \ll k_B T$), both oscillators have equal average thermal energy ($k_B T$) by the equipartition law. If they are nonlinearly coupled, we expect that thermally driven upconversion and downconversion processes are balanced and thermal equilibrium with the reservoir is maintained. However, in the practically relevant nonclassical regime ($\hbar \omega \gtrsim k_B T$), the oscillators have unequal average thermal energies and this raises important fundamental questions. Is the balance between upconversion and downconversion maintained? Are the coupled oscillators in equilibrium with the reservoir? These questions must be reliably answered by the theory before it is used to describe the thermal radiation.

Almost all thermal radiation phenomena (in both classical and nonclassical regimes) are described semi-classically using Kirchhoff’s laws [12,13], fluctuational electrodynamics [14,15] and Langevin coupled mode theory [16,17]. However, we find that none of them can be readily extended for the problem at hand, and the answers are found only within the quantum theory of damping [11,18]. Within the quantum formalism, the problem is non-trivial because of the lack of closed-form analytic solutions and is unsolved in related early works [19–21]. We solve it numerically using the corresponding quantum master equation. Our analysis reveals that the resonators are at thermal equilibrium with the reservoir along with the balance between upconversion and downconversion processes. We also show the failure of alternative semi-classical theory [4,22] which leads to unreasonably large deviation from thermal equilibrium. That comparison further justifies the use of quantum theory for analyzing the nonlinearly coupled passive resonators.

We then use our quantum master equation approach to analyze the thermal emission from the system. We show that by suitably engineering the system to favor the upconversion of thermal radiation from $\omega _1$ to $\omega _2$, one can enhance the thermal emission at $\omega _2$ beyond its linear thermal emission. For a resonant mode at $\omega _2$ of given linewidth, the maximum thermal emission using this nonlinear mechanism is four times larger than the maximum thermal emission achievable in the linear (Planckian) regime. The proposed mechanism realizes the enhancement via spectral redistribution of thermal energy which is fundamentally different from conventional approaches limited to linear media for enhancing far-field [23–25] and near-field [26,27] thermal radiation. We further note that, to guide the design of energy-conversion technologies [2,3], recent interesting works [28–31] have explored bounds on radiative heat transfer in the linear regime. In that context, the finding that the linear bound can be surpassed by the nonlinear upconversion of thermal energy is new and important. Recently, the far-field thermal emission from finite-size subwavelength objects is termed as super-Planckian because it greatly exceeds the flux emitted by a blackbody of the same geometric area. However, since Kirchhoff’s law is still valid and there is no enhancement when the absorption cross-section is used instead of the geometric cross-section, it remains an open question whether this geometric enhancement can be identified as super-Planckian. In the context of this interesting topic [23–25,32], we note that Kirchhoff’s law is not applicable in the present work and the mechanism of nonlinear upconversion paves the way for a regime beyond the super-Planckian enhancement.

Apart from the thermal-emission enhancement, we also find other surprising effects. First, the autocorrelation $g^{(2)}(0)$ of emitted thermal light deviates noticeably from the known value of $g^{(2)}(0)=2$ for thermal light. And second, the emitted thermal radiation exhibits strong correlations between the intensities collected at two distinct frequencies of $\omega _1$ and $\omega _2$. These features are surprising since they are observed by heating the nonlinear system, without using any external signal. From the perspective of quantum optics, we further note that our analysis is markedly different from the analysis of the same system driven by an external coherent pump [19–21]. Most notably, its importance for spectral-engineering thermal radiation [33,34] is not explored. The nascent topic of thermal-radiation engineering using nonlinear media [4–7,35] remains to be explored experimentally. However, we are confident that recent theoretical inquiries including this work will motivate the experiments in the future because of their fundamentally advancing nature and practical utility.

## 2. Quantum theory

We consider a system of two photonic resonant modes of frequencies $\omega _j$ described by operators $\{a_j,a_j^{\dagger }\}$ for $j=[1,2]$ and having frequencies $\omega _1,\omega _2 \approx 2\omega _1$. The resonators contain a $\chi ^{(2)}$ nonlinear medium which facilitates coupling between the resonators via upconversion and downconversion processes. The system Hamiltonian is [19,20]:

The resonators are further coupled linearly with an intrinsic (dissipative) environment/heat bath of harmonic oscillators described by $\{b_{\mathbf {k},d},b_{\mathbf {k},d}^{\dagger }\}$. The radiation from the resonators into an external environment or other channels is modeled by linear coupling with radiation modes which are also harmonic oscillators described by operators $\{b_{\mathbf {k},e},b_{\mathbf {k},e}^{\dagger }\}$. The total Hamiltonian is:

It is straightforward to obtain the temporal dynamics of important relevant quantities from the master equation (Eq. (3)) such as mean photon numbers $\langle a_j^{\dagger } a_j \rangle$ for $j=[1,2]$, where $\langle \cdots \rangle$ denotes the quantum statistical average.

The above expressions also clarify the measurable quantities in an actual experiment. In a practical system, the radiative heat transferred to the external environment at a lower temperature $T_e\;<\;T_d$ can be measured by a suitable detector. By collecting the terms in Eqs. (4) and (5) corresponding to the energy exchange with the external bath, we can quantify the far-field thermal-emission power from the resonators as,

We note that the nonlinear interaction term in the system Hamiltonian given by Eq. (1) commutes with the Hamiltonian of uncoupled oscillators under the frequency matching condition ($\omega _2=2\omega _1$). Therefore, despite the nonlinear coupling, the above solution using the photon number eigenstates is justified. For small mismatch in the frequencies ($\omega _2\approx 2\omega _1$), the solution is approximate. However, we argue that it is reasonable if the frequency mismatch is of the order of nonlinear coupling ($\kappa$). For practical systems, the nonlinearity is very weak ($\kappa \lesssim 10^{-5}\omega _1$ described further below) and for strong nonlinear effects on thermal radiation, the condition $\kappa \sim \gamma$ necessitates the use of resonators of similar linewidths ($\gamma \sim \kappa$). We find that a small frequency mismatch ($|\omega _2-2\omega _1| \sim (\kappa ,\gamma$)) affects the relevant physical quantities with relative deviations of $\mathcal {O}\{\kappa /\omega \} \lesssim 10^{-5}$ which can be considered negligible for practical purpose.

## 3. Thermal equilibrium of resonators

We first analyze the situation where the local resonator heat bath temperature $T_d$ is the same as the temperature of the external environment $T_e$. In the absence of nonlinear coupling $\kappa$, it is straightforward to obtain the steady state density matrix $\rho _{n_1,m_1; n_2,m_2}$ analytically [11,18]. In the photon number basis, it is diagonal ($\rho _{n,m;n,m}$) and is given below:

Figure 2(a) provides an illustration by considering resonant frequencies as $\omega _1=1.57\times 10^{14}$rad/s ($\sim 12\mu$m), $\omega _2=2\omega _1$ ($\sim 6\mu$m) and temperatures $T_d=T_e=T=600$K. The schematic shows the system and various power flows. The mean photon numbers of bath oscillators at these frequencies are $\bar {n}_{\omega _1,T}=0.159$ and $\bar {n}_{\omega _2,T}=0.02$. As demonstrated in this figure, despite the finite nonlinear coupling ($\kappa \neq 0$), the resonator mean photon numbers are equal to those of bath oscillators ($\langle n_{\omega _j} \rangle =\langle a_j^{\dagger } a_j \rangle$). This equality holds irrespective of the temperatures $T$ (higher or lower mean photon numbers of bath oscillators), the coupling parameter $\kappa /\gamma$, and the decay rates $\gamma _{jd}/\gamma$ and $\gamma _{je}/\gamma$. We have normalized these parameters with a suitable $\gamma$ which characterizes the typical linewidth of the resonators. While the nonlinearity is perturbative in comparison with frequencies ($\kappa \ll \omega _j$), we note that the above results are obtained in the strong nonlinear regime ($\kappa \sim \gamma$). As we show in the next section, strong nonlinear effects are indeed observed for the same parameters but under nonequilibrium condition ($T_d \neq T_e$) because of subtle modifications of power flows that otherwise balance each other when the reservoirs are at thermal equilibrium ($T_d = T_e$).

We note that other nonlinear processes of the form $\omega _1+\omega _1'=\omega _2$ where $\omega _1\neq \omega _1'$ are also possible in the presence of $\chi ^{(2)}$ nonlinearity. But, here we focus on extensively studied second-harmonic generation process ($\omega _2=2\omega _1$) from the perspective of potential experiments. If the frequencies are mismatched slightly ($|\omega _2-2\omega _1|\sim \gamma ,\kappa$), the above solution using the photon number eigenstates is approximate as we noted earlier. While very small relative deviations may arise because of the approximation, it can be argued that the oscillators remain at thermal equilibrium with the heat bath. In particular, the nonlinear energy exchange rates given by Eqs. (7) and (8) oscillate sinusoidally over a timescale $\Delta t \sim 1/|\omega _2-2\omega _1|$ in the absence of frequency matching which becomes evident in the interaction picture with respect to the Hamiltonian of uncoupled oscillators. Averaging the flux rates over a sufficiently long duration of time ($\gg \Delta t$), the overall time-averaged energy exchange rates given by Eqs. (7) and (8) are zero. It then follows from Eqs. (4) and (5) that the resonators remain at thermal equilibrium with the bath oscillators. Furthermore, at any given instant of time, there is no violation of energy conservation since all energy flux rates, properly accounted for by the terms in Eqs. (4) and (5), balance each other as shown in the schematic of Fig. 2. We note that a very large frequency mismatch ($|\omega _2-2\omega _1| \gg \gamma ,\kappa$) is beyond the scope of the present theory. It requires consideration of other nonlinear processes such as four-wave mixing to account for the large energy difference given by $\hbar |\omega _2-2\omega _1|$ associated with the microscopic nonlinear mixing process. Also, realizing strong nonlinearity in this regime is highly impractical. Therefore, we do not discuss it here.

As a useful comparison, we further demonstrate that an alternative semiclassical Langevin theory leads to an unreasonable deviation from thermal equilibrium condition for the same parameters. See Fig. 4. We note that semi-classical Langevin theories can accurately describe thermal fluctuations phenomena in the classical regime ($\hbar \omega \ll k_B T$). In the non-classical regime ($\hbar \omega \gtrsim k_B T$), they are reliable for linearly coupled oscillators [16,17] and also for isolated nonlinear (anharmonic) oscillators [4,22]. However, semi-classical Langevin theory fails to describe the nonlinear system considered here. Rather than a simple quantum–classical distinction, this failure stems from its inadequacy in accurately capturing the complicated nonlinear mixing of thermal energy between the oscillators of unequal non-classical frequencies which is required to maintain thermal equilibrium of the oscillators with the heat baths. These important and subtle theoretical details justify the use of quantum theory of damping in the present work which is otherwise uncommon in the field of thermal radiation. The use of quantum theory in this field has been previously considered for studying temporal dynamics of radiative heat transfer for *linearly* coupled plasmonic nanosystems [41]. However, the temporal dynamics can also be studied in the linear regime using the semi-classical Langevin theory [16,17]. Since our focus lies elsewhere, we do not discuss this separate, interesting topic of temporal dynamics here.

## 4. Far-field thermal-emission enhancement

Having established that the quantum theory provides a reasonable description of thermal equilibrium, we now use it to analyze far-field thermal radiation from the same resonant system under the nonequilibrium condition of $T_d \gg T_e$. We assume temperatures $T_d=600$K and $T_e=0$K and compute the thermal radiation power given by Eq. (9). We assume the resonators to be frequency-matched ($\omega _2 =2\omega _1$) but we remark that any frequency mismatch of the order of the linewidth negligibly affects the results described below. Also, throughout the manuscript, the linear coupling between the resonators is considered to be negligible because of the vast difference in the frequencies ($\omega _2 \neq \omega _1, |\omega _2-\omega _1| \gg \gamma _{jl}$).

**Nonlinear upconversion-induced enhancement of thermal emission:** The schematic in Fig. 3(a) depicts various radiative heat exchange channels where the directionalities of the arrows indicate plausible directions of net energy flows under the condition $T_d \gg T_e$. Unlike the thermal equilibrium behavior discussed in the previous section, here the upconversion and downconversion rates are no longer balanced. Under a suitable condition, thermal energy at $\omega _1$ can be upconverted to $\omega _2$ and emitted into the external environment faster than its rate of downconversion back to $\omega _1$. This imbalance provides a new mechanism of enhancing the thermal emission at $\omega _2$. Intuitively, this mechanism is efficient when the spurious decay rates $\gamma _{1e},\gamma _{2d}$ (shown by black arrows in the schematic) are small and favorable decay/coupling rates $\gamma _{1d},\gamma _{2e},\kappa$ (red arrows) are large. Accordingly, we analyze below four possible combinations of these decay rates as shown in the table and use the same color code for all figures in Figs. 3(b), (d), and (e).

Figure 3(b) plots the ratio of thermal emission powers, $F_j = P_j(\kappa )/P_j(0)$ for $j=[1,2]$, as a function of nonlinear coupling $\kappa$. It captures the enhancement of thermal emission beyond the linear regime at $\omega _2$ due to nonlinear upconversion of thermal radiation at $\omega _1$. The thermal emission power in the linear regime is $P_2(0)=2\frac {\gamma _{2e}\gamma _{2d}}{\gamma _2}\hbar \omega _2 \bar {n}_{\omega _2,T_d}$ where $\gamma _2=\gamma _{2e}+\gamma _{2d}$. As evident from the red curve in Fig. 3(b) corresponding to operating regimes favorable for efficient nonlinear upconversion, enhancement factors much greater than $1$ are realized. The additional thermal energy in the emitted power essentially comes from the nonlinear upconversion of thermal energy abundantly available ($\bar {n}_{\omega _1,T} \gg \bar {n}_{\omega _2,T}$) at low frequency $\omega _1$. Consequently, an opposite trend at $\omega _1$ (right figure) is expected and observed which shows the reduction of thermal emission power due to energy depletion.

**Maximum nonlinear enhancements:** We further predict the maximum enhancement compared to linear thermal emission. Since thermal radiation at $\omega _2$ is directly proportional to the mean photon number $\langle n_2 \rangle$, we can maximize it by keeping only the energy exchange channels characterized by decay rate $\gamma _{1d}$ and the nonlinear coupling $\kappa$. Based on this intuition, we numerically find that $\langle n_2 \rangle$ which is otherwise bounded in the linear regime by $\bar {n}_{2}=\gamma _{2d} \bar {n}_{\omega _2,T}/(\gamma _{2d}+\gamma _{2e})$ approaches $\bar {n}_{\omega _2,T}$ in this nonlinear regime. When $\gamma _{2e} \gg \gamma _{2d}$, the linear thermal emission is small since $\bar {n}_{2} \ll \bar {n}_{\omega _2,T}$. For this bad cavity regime, the nonlinear upconversion of thermal energy at $\omega _1$ allows the resonator to reach $\bar {n}_2 = \bar {n}_{\omega _2,T}$, enhancing thermal emission by large factor of $(\gamma _{2e}+\gamma _{2d})/\gamma _{2d}$.

While the above comparison is for given decay rates of the resonant system, we demonstrate another useful comparison in Fig. 3(c) between the maximum thermal emission powers realizable in linear versus nonlinear regimes for a resonant mode of the same total linewidth ($\gamma _2$). It is known that a perfect linear thermal emitter satisfies the rate-matching condition ($\gamma _{2d}=\gamma _{2e}=\gamma _2/2$) [16] and emits the maximum power given by $P_2^{\textrm {max}}(0)= \frac {\gamma _2}{2}\hbar \omega _2 \bar {n}_{\omega _2,T_d}$. This linear limit can be surpassed in the nonlinear regime where the maximum power, $P_2^{\textrm {max}}(\kappa )= 2\gamma _2\hbar \omega _2 \bar {n}_{\omega _2,T_d}$, is obtained when the favorable channels dominate ($\gamma _{1d},\kappa ,\gamma _{2e} \gg \gamma _{2d},\gamma _{1e}$). It follows that the maximum thermal emission obtained via nonlinear upconversion of thermal energy of a *single* resonant mode is four times larger than the maximum thermal emission power in the linear regime. This analysis proves that the known bounds on the enhancement of frequency-selective thermal emission can be surpassed via nonlinear spectral redistribution of energy. This finding is particularly important in the context of recent works on the far-field emission enhancements from subwavelength emitters [23–25] and the fundamental bounds on radiative heat transfer [28–31].

**Modified statistics and biphoton intensity correlations:** In addition to the thermal-emission enhancements beyond the linear regime, we find that the statistical nature of intensity fluctuations quantified by $g^{(2)}(0)$, is also modified because of the nonlinear mixing of thermal energy. Figure 3(d) depicts the modification in $g^{(2)}(0)$ as a function of nonlinear coupling $\kappa$ for both resonators. Evidently, super-bunching $g^{(2)}(0)\;>\;2$ can also be observed for finite nonlinear coupling. Furthermore, Fig. 3(e) shows that the emitted thermal radiation collected by detectors at distinct frequencies $\omega _j$ for $j=[1,2]$ can exhibit correlations. We quantify these correlations as, $C = \frac {\langle I(\omega _1)I(\omega _2) \rangle }{\langle I(\omega _1) \rangle \langle I(\omega _2) \rangle } - 1$, where the correlations are zero when there is no coupling between the resonators ($C=0$ when $\kappa =0$). Here $I(\omega )$ denotes the intensity which is proportional to $\langle a_j^{\dagger }(\omega )a_j(\omega )\rangle$. As demonstrated, the correlations are nonzero for finite nonlinear coupling $\kappa$ and are quite large so as to make them experimentally noticeable. We highlight that these features are obtained by heating the nonlinear system without using any external signal. Also, they indicate new degrees of freedom (statistics and biphoton correlations) for thermal-radiation engineering, which is otherwise concerned with the control of magnitude, directionality, spectrum and polarization of emitted light [33,34]. Finally, we note that we have checked (but not shown) that the nonlinear system does not lead to quantum entanglement (using inseparability criterion) [42] and quantum squeezing [18]. Other interesting questions such as temporal dynamics from the perspective of quantum optics are not explored since the focus is on far-field thermal-radiation which is described in the steady state.

**Potential experimental implementation:** The underlying system considered in the present work is quite well-known [18,43]. Furthermore, its optimization for maximizing the nonlinear frequency mixing is relentlessly explored by many scientists [44–47] because of its importance for biotechnology (bio-sensing) [48], material science (spectroscopy) [49] and quantum science (squeezed, entangled, single photon sources) [18,43]. Our goal here is to point out its potential for thermal science. To accelerate research in that direction, we note experimental observability of these effects by making quantitative predictions based on structures described in the existing literature [44–47]. It was shown in [47] that using typical materials such as AlGaAs ($\chi ^{(2)} \sim 100$pm/V) in optimized microposts and grating structures, the nonlinear coupling $\kappa =\beta \sqrt {\frac {\hbar \omega _1}{2}}$ where $\beta =0.01\frac {\chi ^{(2)}}{\sqrt {\epsilon _0 \lambda _1^{3}}}\frac {2\pi c}{\lambda _1}$ (derivation provided in the appendix) can be realized. For a wavelength of $\lambda _1 = 2\mu$m, this results in the nonlinear coupling $\kappa = 3\times 10^{-8}\omega$. Since dissipative and radiative decay rates should be of the order of the nonlinear coupling $\kappa$ to observe appreciable nonlinear effects, this will require the same resonators of quality factors of $Q \sim 10^{8}$. Another recent work [50] explores the use of multiple-quantum-well semiconductor hetero-structures [51] to realize orders of magnitude larger material nonlinearities ($\chi ^{(2)} \sim 10^{5}$pm/V). A combination of these two approaches can potentially realize $\kappa \sim 10^{-5}\omega$ (as explored in Fig. 2 and Fig. 3), lowering the quality factor requirements to $Q \sim 10^{5}$ to observe the predicted nonlinear effects. One possibility to observe strong nonlinear effects with lower quality factor resonators is to harness the large density of surface polaritonic states ($Q \gtrsim 100$) of planar media as explored in our earlier work [6]. An extension of the present work for that geometry requires additional considerations because of many resonant modes and therefore, it will be analyzed in our subsequent work. We finally summarize all suitable alternatives for detecting nonlinear thermal radiation effects in passive systems. For the second-order nonlinear process considered here, this includes large density of highly confined phonon polaritons in the near-field of polaritonic media [6,52], multiple-quantum-well hetero-structures for giant nonlinearities ($\chi ^{(2)} \sim 10^{5}$pm/V) [50], enhanced nonlinearities in two-dimensional materials such as bilayer graphene (tunable $\chi ^{(2)} \sim 10^{5}$pm/V) [53], high-Q bound states in the continuum [54,55] and inverse-design optimization [56]. We are confident that an optimized nonlinear system designed in the future will reveal the predicted effects and also be useful for many other applications.

## 5. Conclusion

We analyzed thermal radiation from a system of two *nonlinearly* coupled resonators of *distinct* frequencies which also requires careful examination of its thermal equilibrium behavior. At frequencies ($\hbar \omega \gtrsim k_B T$), the oscillators have unequal average thermal energies and it is not obvious whether the coupled oscillators are at thermal equilibrium with the reservoir. Our quantum theoretic approach describes the thermal equilibrium behavior of the coupled oscillators reasonably well. We also note that it cannot be adequately captured by semi-classical Langevin theory which is otherwise useful in the linear regime [16,17] and for isolated nonlinear (anharmonic) oscillators [4,22]. These theoretical findings are generalizable to other systems described using such nonlinear harmonic oscillators and invite similar studies of nonlinear thermal-fluctuations phenomena in other fields of research e.g. thermal nonlinearities in passive (undriven) mechanical oscillators [57].

Most importantly, we provided a new mechanism of enhancing the far-field thermal emission beyond the linear blackbody limits via nonlinear upconversion. Spectral-engineering thermal radiation using nonlinear media is a nascent topic and the finding that it can allow one to surpass the linear Planckian bounds is fundamentally important in view of recent inquiries concerning super-Planckian thermal emission from subwavelength emitters [23–25,32] and bounds on radiative heat transfer [28–31]. While we did not establish any tight bounds on feasible nonlinear enhancements, our preliminary results motivate future work on thermal photonic nonlinearities in other systems and inquiries of generalized Kirchhoff’s laws [58] and thermal-emission sum rules for nonlinear media. We also found that because of the nonlinear mixing, the emitted thermal light exhibits nontrivial statistics ($g^{(2)}(0) \neq 2$) and biphoton intensity correlations. These features are surprising from the perspective of quantum optics because they can be observed by heating a nonlinear system without the need for any external signal. Finally, from the perspective of theory development in the field of thermal radiation, we note that recent scientific efforts have extended the existing fluctuational electrodynamic paradigm to explore non-traditional nonreciprocal [59,60] and nonequilibrium [32,61,62] systems. Our present work aspires to go beyond these regimes to study nonlinearities and while doing so, also departs from the traditional, semi-classical theoretical paradigms. We are confident that additional fundamental discoveries will be made in this largely unexplored territory in the near future.

## Appendix

We analyze the same doubly resonant system containing a $\chi ^{(2)}$ nonlinear medium and supporting resonances at frequencies $\omega _1$ and $\omega _2=2\omega _1$. We use the semi-classical Langevin theory and show that it leads to large deviation from thermal equilibrium condition, which is otherwise captured reasonably well using the quantum theory of damping described in the main text.

**Semi-classical theory:** We directly introduce the Langevin (temporal coupled mode) equations for the resonator amplitudes based on previous work [4]. By denoting the field amplitude as $a_j(t)$ for $j=[1,2]$ normalized such that $|a_j|^{2}$ denotes the mode energy, we write the following coupled mode equations:

We now obtain the FD coefficients $D_j$ by transforming the stochastic ODE into corresponding Fokker-Planck equation for the probability distribution $P(a_1,a_1^{*},a_2,a_2^{*})$ which is:

We numerically simulate the stochastic equations in Eqs. (12) and (13) and compute the mean photon numbers in the steady state as $\langle n_{\omega _j} \rangle = \frac {\langle |a_j(t)|^{2} \rangle }{\hbar \omega _j}$. Figure 4 compares the semi-classical theory with quantum theory for the exact same parameters as considered in Fig. 2(a) of the main text. In the figure, we have used the relation given by Eq. (18) to demonstrate the dependence of the mean photon numbers on the nonlinear coupling $\kappa$ instead of the nonlinear coefficient $\beta$. Evidently, the steady state values obtained using the semi-classical theory indicate significant deviation from equilibrium mean photon numbers of heat-bath oscillators as a function of the nonlinear coupling. Since a large relative deviation from thermal equilibrium occurs despite perturbative nature of the nonlinearity ($\kappa \approx 10^{-5}\omega$) and the perfect frequency matching condition $\omega _2=2\omega _1$, we conclude that the Langevin theory fails to accurately describe this particular nonlinear system. Intuitively, this may be expected since the resonant frequencies are chosen in the non-classical regime ($\hbar \omega \gtrsim k_B T$). However, it is not obvious because the semi-classical theory is otherwise reliable in this regime for linearly coupled resonators [16,17] and for isolated single nonlinear (anharmonic) oscillators [4,22]. We therefore conclude that the semi-classical theory fails because it cannot accurately capture the complicated nonlinear mixing of thermal energy between the resonant modes of non-classical frequencies which can ensure thermal equilibrium of resonant modes with the heat baths. Similar failures of Langevin theory applied to nonlinear systems (as opposed to quantum versus classical distinctions) have been noted before by van Kampen [64].

We note that we have also explored several modified forms of the semi-classical Langevin theory unsuccessfully after our prior work [4] until now, to describe the thermal equilibrium behavior of nonlinearly coupled resonators. One approach involves considering multiplicative noise where the terms $D_j$ in the coupled mode equations [Eqs. (12) and (13)] are dependent on mode amplitudes $a_j$ for $j=[1,2]$. Their unknown functional form should be determined using the Fokker-Planck equation by enforcing the known equilibrium probability distribution. However, the complicated form of the multiplicative noise solution lacks any physical justification [64] and also does not yield a consistent theory. Nonetheless, we leave it as an open question whether thermal equilibrium behavior of the nonlinear system considered here can be reproduced reliably using a semiclassical theory.

## Funding

Defense Advanced Research Projects Agency (N66001-17-1-4048, HR00111820046); Purdue University (Lillian Gilbreth Postdoctoral Fellowship Program); National Science Foundation (DMR-1454836); Cornell Center for Materials Research (DMR1719875).

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **D. Chang, V. Vuletić, and M. Lukin, “Quantum nonlinear optics–photon by photon,” Nat. Photonics **8**(9), 685–694 (2014). [CrossRef]

**2. **S. Fan, “Thermal photonics and energy applications,” Joule **1**(2), 264–273 (2017). [CrossRef]

**3. **E. Tervo, E. Bagherisereshki, and Z. Zhang, “Near-field radiative thermoelectric energy converters: a review,” Front. Energy **12**(1), 5–21 (2018). [CrossRef]

**4. **C. Khandekar, A. Pick, S. Johnson, and A. Rodriguez, “Radiative heat transfer in nonlinear kerr media,” Phys. Rev. B **91**(11), 115406 (2015). [CrossRef]

**5. **C. Khandekar, Z. Lin, and A. Rodriguez, “Thermal radiation from optically driven kerr (*χ* (3)) photonic cavities,” Appl. Phys. Lett. **106**(15), 151109 (2015). [CrossRef]

**6. **C. Khandekar and A. Rodriguez, “Near-field thermal upconversion and energy transfer through a kerr medium,” Opt. Express **25**(19), 23164–23180 (2017). [CrossRef]

**7. **C. Khandekar, R. Messina, and A. Rodriguez, “Near-field refrigeration and tunable heat exchange through four-wave mixing,” AIP Adv. **8**(5), 055029 (2018). [CrossRef]

**8. **R. W. Boyd, * Nonlinear optics* (Elsevier, 2003).

**9. **S. E. Harris, J. E. Field, and A. Imamoğlu, “Nonlinear optical processes using electromagnetically induced transparency,” Phys. Rev. Lett. **64**(10), 1107–1110 (1990). [CrossRef]

**10. **H. Schmidt and A. Imamoglu, “Giant kerr nonlinearities obtained by electromagnetically induced transparency,” Opt. Lett. **21**(23), 1936–1938 (1996). [CrossRef]

**11. **H. Breuer and F. Petruccione, * The theory of open quantum systems* (Oxford University, 2002).

**12. **L. Landau and E. Lifshitz, “Course of theoretical physics, volume 5,” Publ. Butterworth-Heinemann **3** (1980).

**13. **C. Luo, A. Narayanaswamy, G. Chen, and J. Joannopoulos, “Thermal radiation from photonic crystals: a direct calculation,” Phys. Rev. Lett. **93**(21), 213905 (2004). [CrossRef]

**14. **S. Rytov, “Theory of electric fluctuations and thermal radiation,” AFCRC-TR **59**, 162 (1959).

**15. **C. Otey, L. Zhu, S. Sandhu, and S. Fan, “Fluctuational electrodynamics calculations of near-field heat transfer in non-planar geometries: A brief overview,” J. Quant. Spectrosc. Radiat. Transfer **132**, 3–11 (2014). [CrossRef]

**16. **L. Zhu, S. Sandhu, C. Otey, S. Fan, M. Sinclair, and T. Shan Luk, “Temporal coupled mode theory for thermal emission from a single thermal emitter supporting either a single mode or an orthogonal set of modes,” Appl. Phys. Lett. **102**(10), 103104 (2013). [CrossRef]

**17. **A. Karalis and J. Joannopoulos, “Temporal coupled-mode theory model for resonant near-field thermophotovoltaics,” Appl. Phys. Lett. **107**(14), 141108 (2015). [CrossRef]

**18. **M. Scully and M. Zubairy, * Quantum optics* (AAPT, 1999).

**19. **G. Agarwal, “Quantum theory of second harmonic generation,” Opt. Commun. **1**(3), 132–134 (1969). [CrossRef]

**20. **P. Drummond and M. Hillery, * The quantum theory of nonlinear optics* (Cambridge University, 2014).

**21. **M. Kozierowski and R. Tanaś, “Quantum fluctuations in second-harmonic light generation,” Opt. Commun. **21**(2), 229–231 (1977). [CrossRef]

**22. **M. Dykman and M. Krivoglaz, “Spectral distribution of nonlinear oscillators with nonlinear friction due to a medium,” Phys. Status Solidi B **68**(1), 111–123 (1975). [CrossRef]

**23. **S.-A. Biehs and P. Ben-Abdallah, “Revisiting super-planckian thermal emission in the far-field regime,” Phys. Rev. B **93**(16), 165405 (2016). [CrossRef]

**24. **V. Fernández-Hurtado, A. Fernández-Domínguez, J. Feist, F. García-Vidal, and J. Cuevas, “Super-planckian far-field radiative heat transfer,” Phys. Rev. B **97**(4), 045408 (2018). [CrossRef]

**25. **D. Thompson, L. Zhu, R. Mittapally, S. Sadat, Z. Xing, P. McArdle, M. M. Qazilbash, P. Reddy, and E. Meyhofer, “Hundred-fold enhancement in far-field radiative heat transfer over the blackbody limit,” Nature **561**(7722), 216–221 (2018). [CrossRef]

**26. **J. Yang, W. Du, Y. Su, Y. Fu, S. Gong, S. He, and Y. Ma, “Observing of the super-planckian near-field thermal radiation between graphene sheets,” Nat. Commun. **9**(1), 4033 (2018). [CrossRef]

**27. **K. Kim, B. Song, V. Fernández-Hurtado, W. Lee, W. Jeong, L. Cui, D. Thompson, J. Feist, M. Reid, F. García-Vidal, E. Meyhofer, and P. Reddy, “Radiative heat transfer in the extreme near field,” Nature **528**(7582), 387–391 (2015). [CrossRef]

**28. **O. D. Miller, A. G. Polimeridis, M. H. Reid, C. W. Hsu, B. G. DeLacy, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, “Fundamental limits to optical response in absorptive systems,” Opt. Express **24**(4), 3329–3364 (2016). [CrossRef]

**29. **S. Buddhiraju, P. Santhanam, and S. Fan, “Thermodynamic limits of energy harvesting from outgoing thermal radiation,” Proc. Natl. Acad. Sci. **115**(16), E3609–E3615 (2018). [CrossRef]

**30. **S. Molesky, W. Jin, P. Venkataram, and A. Rodriguez, “Bounds on absorption and thermal radiation for arbitrary objects,” arXiv:1907.04418 (2019).

**31. **S. Molesky, P. Venkataram, W. Jin, and A. Rodriguez, “Fundamental limits to radiative heat transfer: theory,” arXiv:1907.03000 (2019).

**32. **J.-J. Greffet, P. Bouchon, G. Brucoli, and F. Marquier, “Light emission by nonequilibrium bodies: local kirchhoff law,” Phys. Rev. X **8**(2), 021008 (2018). [CrossRef]

**33. **W. Li and S. Fan, “Nanophotonic control of thermal radiation for energy applications,” Opt. Express **26**(12), 15995–16021 (2018). [CrossRef]

**34. **D. Baranov, Y. Xiao, I. Nechepurenko, A. Krasnok, A. Alù, and M. Kats, “Nanophotonic engineering of far-field thermal emitters,” Nat. Mater. **18**(9), 920–930 (2019). [CrossRef]

**35. **H. Soo and M. Krüger, “Fluctuational electrodynamics for nonlinear materials in and out of thermal equilibrium,” Phys. Rev. B **97**(4), 045412 (2018). [CrossRef]

**36. **J. Joannopoulos, S. Johnson, J. Winn, and R. Meade, * Photonic crystals: Molding the flow of light* (Princeton University, 2011).

**37. **J. Bravo-Abad, S. Fan, S. Johnson, J. Joannopoulos, and M. Soljacic, “Modeling nonlinear optical phenomena in nanophotonics,” J. Lightwave Technol. **25**(9), 2539–2546 (2007). [CrossRef]

**38. **A. Rodriguez, M. Soljačić, J. Joannopoulos, and S. Johnson, “*χ* (2) and *χ* (3) harmonic generation at a critical power in inhomogeneous doubly resonant cavities,” Opt. Express **15**(12), 7303–7318 (2007). [CrossRef]

**39. **L.-P. Yang and Z. Jacob, “Engineering first-order quantum phase transitions for weak signal detection,” arXiv preprint arXiv:1905.07420 (2019).

**40. **L.-P. Yang and Z. Jacob, “Quantum critical detector: amplifying weak signals using discontinuous quantum phase transitions,” Opt. Express **27**(8), 10482–10494 (2019). [CrossRef]

**41. **S.-A. Biehs and G. Agarwal, “Dynamical quantum theory of heat transfer between plasmonic nanosystems,” J. Opt. Soc. Am. B **30**(3), 700–707 (2013). [CrossRef]

**42. **A. Peres, “Separability criterion for density matrices,” Phys. Rev. Lett. **77**(8), 1413–1415 (1996). [CrossRef]

**43. **D. Walls and G. Milburn, * Quantum optics* (Springer, 2008).

**44. **W. Pernice, C. Xiong, C. Schuck, and H. Tang, “Second harmonic generation in phase matched aluminum nitride waveguides and micro-ring resonators,” Appl. Phys. Lett. **100**(22), 223501 (2012). [CrossRef]

**45. **K. Rivoire, Z. Lin, F. Hatami, W. Masselink, and J. Vučković, “Second harmonic generation in gallium phosphide photonic crystal nanocavities with ultralow continuous wave pump power,” Opt. Express **17**(25), 22609–22615 (2009). [CrossRef]

**46. **Z.-F. Bi, A. Rodriguez, H. Hashemi, D. Duchesne, M. Loncar, K.-M. Wang, and S. Johnson, “High-efficiency second-harmonic generation in doubly-resonant *χ* (2) microring resonators,” Opt. Express **20**(7), 7526–7543 (2012). [CrossRef]

**47. **Z. Lin, X. Liang, M. Lončar, S. Johnson, and A. Rodriguez, “Cavity-enhanced second-harmonic generation via nonlinear-overlap optimization,” Optica **3**(3), 233–238 (2016). [CrossRef]

**48. **P. Campagnola and L. Loew, “Second-harmonic imaging microscopy for visualizing biomolecular arrays in cells, tissues and organisms,” Nat. Biotechnol. **21**(11), 1356–1360 (2003). [CrossRef]

**49. **T. Heinz, C. Chen, D. Ricard, and Y. Shen, “Spectroscopy of molecular monolayers by resonant second-harmonic generation,” Phys. Rev. Lett. **48**(7), 478–481 (1982). [CrossRef]

**50. **J. Lee, M. Tymchenko, C. Argyropoulos, P.-Y. Chen, F. Lu, F. Demmerle, G. Boehm, M.-C. Amann, A. Alu, and M. Belkin, “Giant nonlinear response from plasmonic metasurfaces coupled to intersubband transitions,” Nature **511**(7507), 65–69 (2014). [CrossRef]

**51. **E. Rosencher, A. Fiore, B. Vinter, V. Berger, P. Bois, and J. Nagle, “Quantum engineering of optical nonlinearities,” Science **271**(5246), 168–173 (1996). [CrossRef]

**52. **N. Rivera, G. Rosolen, J. Joannopoulos, I. Kaminer, and M. Soljačić, “Making two-photon processes dominate one-photon processes using mid-ir phonon polaritons,” Proc. Natl. Acad. Sci. **114**(52), 13607–13612 (2017). [CrossRef]

**53. **S. Wu, L. Mao, A. Jones, W. Yao, C. Zhang, and X. Xu, “Quantum-enhanced tunable second-order optical nonlinearity in bilayer graphene,” Nano Lett. **12**(4), 2032–2036 (2012). [CrossRef]

**54. **L. Carletti, K. Koshelev, C. De Angelis, and Y. Kivshar, “Giant nonlinear response at the nanoscale driven by bound states in the continuum,” Phys. Rev. Lett. **121**(3), 033903 (2018). [CrossRef]

**55. **M. Minkov, D. Gerace, and S. Fan, “Doubly resonant *χ*^{(2)} nonlinear photonic crystal cavity based on a bound state in the continuum,” Optica **6**(8), 1039–1045 (2019). [CrossRef]

**56. **C. Sitawarin, W. Jin, Z. Lin, and A. Rodriguez, “Inverse-designed photonic fibers and metasurfaces for nonlinear frequency conversion,” Photonics Res. **6**(5), B82–B89 (2018). [CrossRef]

**57. **J. Gieseler, L. Novotny, and R. Quidant, “Thermal nonlinearities in a nanomechanical oscillator,” Nat. Phys. **9**(12), 806–810 (2013). [CrossRef]

**58. **D. Miller, L. Zhu, and S. Fan, “Universal modal radiation laws for all thermal emitters,” Proc. Natl. Acad. Sci. **114**(17), 4336–4341 (2017). [CrossRef]

**59. **L. Zhu, Y. Guo, and S. Fan, “Theory of many-body radiative heat transfer without the constraint of reciprocity,” Phys. Rev. B **97**(9), 094302 (2018). [CrossRef]

**60. **C. Khandekar and Z. Jacob, “Thermal spin photonics in the near-field of nonreciprocal media,” New J. Phys. **21**(10), 103030 (2019). [CrossRef]

**61. **W. Jin, A. Polimeridis, and A. Rodriguez, “Temperature control of thermal radiation from composite bodies,” Phys. Rev. B **93**(12), 121403 (2016). [CrossRef]

**62. **C. Khandekar and Z. Jacob, “Circularly polarized thermal radiation from nonequilibrium coupled antennas,” Phys. Rev. Appl. **12**(1), 014053 (2019). [CrossRef]

**63. **F. Moss and P. McClintock, * Noise in nonlinear dynamical systems*, vol. 2 (Cambridge University, 1989).

**64. **N. Van Kampen, * Stochastic processes in physics and chemistry*, vol. 1 (Elsevier, 1992).