Quantum nonlinear mixing of thermal photons to surpass the blackbody limit

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


I. 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 χ (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 χ (2) nonlinear medium and discover new fundamental possibilities for the field of thermal radiation.
Our system depicted in figure 1 consists of two resonators at frequencies ω 1 and ω 2 ≈ 2ω 1 , and contains a χ (2) nonlinear material. In stark contrast to prior works in the quantum-optics literature where an external drive is used [8][9][10], here the resonators are driven by low power thermal fluctuations inside the medium. Because of χ (2) nonlinearity, two photons at ω 1 can get upconverted to yield a single photon at ω 2 , and a single photon at ω 2 can get downconverted to yield two photons at ω 1 . The resonantly enhanced, thermally driven upconversion and downconversion processes can alter the well-known Planck's distribution as illustrated in figure 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 Θ(ω, T ) = ω/[exp( ω/k B T ) − 1], above the vacuum zero point energy [11]. In the classical regime ( ω 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 ( ω 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.
However, we find that none of them can be readily extended for the problem at hand, and the answers are found only within ii ω 2 ≈ 2ω 1 1 10 100 χ (2) Upconversion Wavelength (μm)  (2) Reshaping blackbody distribution through nonlinear mixing of thermal photons FIG. 1. We rigorously analyze thermal radiation from a resonant system depicted in the inset. It is comprised of two resonators of frequencies ω1 and ω2 ≈ 2ω1, and contains a χ (2) nonlinear medium. The resonantly enhanced nonlinear interaction will modify the Planck's blackbody distribution, potentially allowing enhancement of thermal emission at wavelengths where it is otherwise exponentially suppressed.
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][20][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, when a certain frequency matching condition (ω 2 = 2ω 1 ) is strictly satisfied. For imperfectly matched resonators (ω 2 = 2ω 1 ), a small yet reasonable perturbation from thermal equilibrium condition is observed which follows from the microscopic nonlinear process. We also show the failure of alternative semi-classical theory [4,22] which leads to unreasonably large deviation from thermal equilibrium despite the frequency matching condition. 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 ω 1 to ω 2 , one can enhance the thermal emission at ω 2 beyond its linear thermal emission. For a resonant mode at ω 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][24][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][29][30][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][24][25]32], we note that Kirchhoff's law is not applicable in present work and the mechanism of nonlinear upconversion paves the way for a regime beyond 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 ω 1 and ω 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][20][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][5][6][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.

II. QUANTUM THEORY
We consider a system of two photonic resonant modes of frequencies ω j described by operators {a j , a † j } for j = [1, 2] and having frequencies ω 1 , ω 2 ≈ 2ω 1 . The resonators contain a χ (2) nonlinear medium which facilitates coupling between the resonators via upconversion and downconversion processes. The system Hamiltonian is [19,20]: where κ denotes weak nonlinear coupling (κ ω) between the modes. The connection of the above phenomenological Hamiltonian with Maxwell's equations can be made by comparing with corresponding classical coupled mode theory widely used in nanophotonics and derived using Maxwell's equations [36]. The resonant frequencies are calculated by solving the eigenmodes of Maxwell's equations while decay and other coupling iii rates are obtained using perturbation theory [37,38]. Using this approach, we show in the supplement that the nonlinear coupling κ depends on the overlap integral of linear resonant mode profiles.
The resonators are further coupled linearly with an intrinsic (dissipative) environment/heat bath of harmonic oscillators described by {b k,d , b † k,d }. 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 k,e , b † k,e }. The total Hamiltonian is: [1,2] Using the standard techniques of analyzing open quantum systems [11,18], we solve the reduced dynamics of the two oscillators by tracing over the heat bath degrees of freedom. In the regime of weak coupling g k,l,j with a continuum of bath oscillators, we assume that the initial state is ρ ⊗ ρ B where ρ is a joint density operator of coupled resonators and ρ B is an equilibrium density operator of bath oscillators. Using Markov approximation and Gibbs distribution for bath oscillators, we obtain the following simplified master equation in the Schrodinger picture: [1,2] γ j,l (n ωj , The above well-known form captures the interaction of each resonator with heat bath/environment using the associated decay rates γ j,l . Here, j = [1, 2] denotes the resonator mode and l = [d, e] denotes either dissipative (d) heat bath or extrinsic (e) radiative environment at temperature T l . For a practical system, the decay rate is related to the quality factor Q of the resonator via the relation γ = ω/2Q. Also,n ω, is the mean photon number of a harmonic oscillator of frequency ω at thermodynamic equilibrium temperature T . It is straightforward to obtain the temporal dynamics of important relevant quantities from the master equation (3) such as mean photon numbers a † j a j for j = [1,2], where · · · denotes the quantum statistical average.
From these equations, one can obtain the rates of energy loss for each mode by multiplying (4) by ω 1 and (5) by ω 2 . In the steady state, it then follows that the total power flow from mode j = [1, 2] to bath l = [d, e] is: Similarly, the total power flows leaving the modes via nonlinear coupling κ are: The equality P 2→1 = −P 1→2 when ω 2 = 2ω 1 indicates the energy conservation satisfied by the underlying microscopic nonlinear process. Deviation from this frequency matching condition has important implications for equilibrium mean photon numbers of the oscillators as we explain further below. 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. We can thus quantify the far-field thermal emission power from the resonators as, for j = [1,2]. The calculation of these relevant quantities in the nonlinear regime (κ = 0) is however not as straightforward as it is for linear systems. Although quantum Langevin equations for the operators can be written in the nonlinear regime, they cannot be analytically solved. Also, the equations of motion for operators such as a † 2 1 a 2 obtained from the master equation contain higher order terms such as a † 3 1 a 1 a 2 , leading to an infinite set of equations where higher order terms are not necessarily negligible. Therefore, we directly obtain the steady state density matrix from the master equation by numerically solving it [39,40]. We write the master equation (3) in its matrix form in the basis of photon number states. The vector form of the joint density operator is (see Appendix B in [39]) iv with |n 1 , m 1 ; n 2 , m 2 = (|n 1 ⊗ n 2 |) ⊗ (|m 1 ⊗ m 2 |). n denotes number of photons at ω 1 and m denotes number of photons at ω 2 . It follows from the trace condition (Trρ = 1) that n m ρ n,m;n,m = 1. We obtain the steady state solution of the density matrix by solving the master equation (3) and the trace condition together. We work in the regime where this problem can be solved numerically by introducing a cutoff n c limiting the upper bounds in (10) such that all the terms ρ n1,m1;n2,m2 beyond the cut-off (n j , m j > n c ) are negligible and can be safely ignored.

III. 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 κ, it is straightforward to obtain the steady state density matrix ρ n1,m1;n2,m2 analytically [11,18]. In the photon number basis, it is diagonal (ρ n,m;n,m ) and is given below: In the presence of nonlinear coupling κ, by employing our numerical approach, we find that the same steady state solution (with zero off-diagonal elements) satisfies the master equation (3) when the frequency matching condition ω 2 = 2ω 1 is strictly satisfied. The linear solution leads to zero nonlinear correction to the energy of the resonant system, Tr[ρ(κa † 2 1 a 2 + κ * a 2 1 a † 2 )] = 0, despite the nonlinear coupling (κ = 0). Consequently, the nonlinear terms in (4) and (5) vanish leading to thermal equilibrium with bath oscillators in the strictest sense ( a † j a j =n ωj ,T ). However, we caution that for imperfectly matched resonators (ω 2 = 2ω 1 ), the density matrix is no longer diagonal and the nonlinear correction to the energy is nonzero due to off-diagonal entries in the density matrix. We carefully analyze this situation further below. Figure 2 provides an illustration by considering resonant frequencies as ω 1 = 1.57 × 10 14 rad/s (∼ 12µm), ω 2 = 2ω 1 (∼ 6µm) and temperatures T d = T e = T = 600K. The inset of figure 2(a) schematically shows the system and various power flows. The mean photon numbers of bath oscillators at these frequencies aren ω1,T = 0.159 andn ω2,T = 0.02. Figure 2(a) demonstrates that despite the finite nonlinear coupling (κ = 0), the resonator mean photon numbers are equal to those of bath oscillators ( n ωj = a † j a j ). This equality holds irrespective of the temperatures T (higher or lower mean photon numbers of bath oscillators), the coupling parameter κ/γ, and the decay rates γ jd /γ and γ je /γ. We have normalized these parameters with a suitable γ which  2. (a) Steady state mean photon numbers ( nω j ) of nonlinearly coupled resonators (κ = 0) are equal to equilibrium mean photon numbers (nω j ,T ) of bath oscillators (T d = Te = T ) irrespective of nonlinear coupling and decay rates. The equality holds when the frequency matching condition, ω2 = 2ω1, is strictly satisfied. The steady state density matrix is diagonal and same as in the linear regime given by Eq.11 (b) In presence of detuning (ω2 − 2ω1 ∼ γ), the density matrix has nonzero off-diagonal elements (top left). However, their contribution is orders of magnitude small compared to diagonal elements (top right). Nonetheless, this causes small relative deviations (∆n/n ∼ 10 −5 ) from equilibrium mean photon numbers of both resonators as illustrated by bottom two panels. These deviations and their dependence on detuning follow from the underlying nonlinear process as explained in the main text. For illustration, we have assumed all decay rates to be equal to γ/ω1 = 10 −5 . We emphasize that the results are obtained for strong nonlinearity (κ ∼ γ). Although nonlinear effects are apparently negligible when T d = Te as considered here, we demonstrate in fig.3 that the same parameters lead to strong nonlinear effects under nonequilibrium condition T d Te because of the subtle modifications of power flows.
characterizes the typical linewidth of the resonators. While the nonlinearity is perturbative in comparison with frequencies (κ ω j ), we note that the above results are obtained in the strong nonlinear regime (κ ∼ γ). As we show in the next section, strong nonlinear effects are indeed observed for the same parameters but under v nonequilibrium condition (T d = 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 now rigorously analyze thermal equilibrium behavior of resonators when the frequencies are mismatched by small detuning (ω 2 − 2ω 1 ∼ γ j ), as promised. Figure 2(b) examines such a system for realistic nonlinear coupling κ/ω ∼ 10 −5 . To observe strong nonlinear effects, small linewidths γ ∼ κ are required.
Therefore, and also for the purpose of illustration, we assume all decay rates to be equal to decay rate γ = 10 −5 ω 1 and also normalize the nonlinear coupling (κ/γ). The top left figure demonstrates the sparsity of the numerically computed density matrix solution. As illustrated, certain off-diagonal elements are nonzero for imperfectly matched resonators. However, the contribution of these off-diagonal terms is orders of magnitude smaller than the dominant diagonal terms. The top right figure in 2(b) compares the density matrix elements ρ 0,0;0,0 , ρ 0,1;0,1 , ρ 1,0;1,0 and ρ 0,1;2,0 as a function of detuning |ω 2 − 2ω 1 |/γ on a logarithmic scale. Evidently, the nonzero off-diagonal element ρ 0,1;2,0 is more than five orders of magnitude smaller than the dominant diagonal terms.
We find that the frequency mismatch leads to small deviation from thermal equilibrium with bath oscillators as shown by the bottom two figures in fig.2(b). The figures display the relative deviation from equilibrium mean photon numbers when (ω 2 − 2ω 1 ) ∼ γ. The deviation arises because of the non-energy-conserving microscopic nonlinear process. When ω 2 = 2ω 1 , the power lost via nonlinear coupling by mode a 1 given by (7) is not equal to the power gained by mode a 2 given by (8) and this difference is compensated by heat exchange with respective baths leading to deviation from thermal equilibrium. For instance, when 2ω 1 > ω 2 , it follows from the underlying microscopic process that the power provided by mode a 1 is larger than that accepted by mode a 2 . The redundant power then goes to the heat bath which implies that a † 1 a 1 >n ωj ,T or ∆n ω1 > 0 as observed in the bottom left figure of fig 2(b). Similar microscopic understanding also explains the opposite deviation at ω 2 shown in the right figure and when 2ω 1 < ω 2 shown in both figures. A large detuning of the order of frequency (|ω 2 −2ω 1 | ∼ ω j γ j ) can lead to large deviation from thermal equilibrium. However, that situation is beyond the scope of the present theory since the accurate description of the highly non-energy conserving nonlinear process (ω 2 ≈ 2ω 1 ) requires additional considerations. Also, realizing strong nonlinearity (κ ∼ γ) in this regime is impractical and therefore, we do not discuss it here. We note that similar deviation is observed for two linearly coupled harmonic oscillators of distinct frequencies. However, for linearly coupled oscillators, one can transform them into normal modes which follow the same harmonic-oscillators commutation algebra and realize thermal equilibrium of new normal modes in the strictest sense.
For nonlinearly coupled oscillators, such a mathematical transformation of operators does not exist. Nonetheless, since the nonlinearity is perturbative and small in practical systems (κ ω j ), the above description which allows and explains small relative deviations (∆n/n ∼ 10 −5 ) of the order of the perturbation (κ ≈ 10 −5 ω) is adequate.
In the supplement, we show that an alternative semiclassical Langevin theory leads to an unreasonable deviation from thermal equilibrium for the same parameters despite the frequency matching condition (ω 2 = 2ω 1 ). We note that semi-classical Langevin theories can accurately describe thermal fluctuations phenomena in the classical regime ( ω k B T ). In the non-classical regime ( ω k B T ), they are reliable for linearly coupled oscillators [16,17] and also for isolated nonlinear (anharmonic) oscillators [4,22]. But, they cannot be used for nonlinearly coupled oscillators of distinct frequencies. These important 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 here.

IV. 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 T e . We assume temperatures T d = 600K and T e = 0K and compute the thermal radiation power given by (9). We assume the resonators to be frequency-matched (ω 2 = 2ω 1 ) but note that any frequency mismatch of the order of the linewidth negligibly affects the results described below. Our emphasis on imperfect frequency matching in the previous section was primarily meant for rigorous analysis of thermal equilibrium and is not repeated here. Also, throughout the manuscript, the linear coupling between the resonators is considered to be negligible because of the vast difference in the frequencies (ω 2 = ω 1 , |ω 2 − ω 1 | γ jl ). Nonlinear upconversion-induced enhancement of thermal emission: The schematic in figure 3(a) vi depicts various radiative heat exchange channels where the directionalities of the arrows indicate plausible directions of net energy flows under the condition T d 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 ω 1 can be upconverted to ω 2 and emitted into the external environment faster than its rate of downconversion back to ω 1 . This imbalance provides a new mechanism of enhancing the thermal emission at ω 2 . Intuitively, this mechanism is efficient when the spurious decay rates γ 1e , γ 2d (shown by black arrows in the schematic) are small and favorable decay/coupling rates γ 1d , γ 2e , κ (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 3(b,c,d). Figure 3(b) plots the ratio of thermal emission powers, F j = P j (κ)/P j (0) for j = [1,2], as a function of nonlinear coupling κ. It captures the enhancement of thermal emission beyond the linear regime at ω 2 due to nonlinear upconversion of thermal radiation at ω 1 . The thermal emission power in the linear regime is P 2 (0) = 2 γ2eγ 2d γ2 ω 2nω2,T d where γ 2 = γ 2e + γ 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 (n ω1,T n ω2,T ) at low frequency ω 1 . Consequently, an opposite trend at ω 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 ω 2 is directly proportional to the mean photon number n 2 , we can maximize it by keeping only the energy exchange channels characterized by decay rate γ 1d and the nonlinear coupling κ. Based on this intuition, we numerically find that n 2 which is otherwise bounded in the linear regime byn 2 = γ 2dnω2,T /(γ 2d + γ 2e ) approachesn ω2,T in this nonlinear regime.
While the above comparison is for given decay rates of the resonant system, we demonstrate another useful comparison in the inset of fig.3(b) between the maximum thermal emission powers realizable in linear versus nonlinear regimes for a resonant mode of the same total linewidth (γ 2 ). It is known that a perfect linear (a) Far-field thermal emission from a doubly resonant nanophotonic system at temperature T d = 600K into external environment at temperature Te = 0K. Various decay/coupling channels are shown as arrows with their directionalities indicating the flow of energy in the regime T d Te.
Red arrows indicate favorable channels for nonlinear-upconversion-induced enhancement of thermal emission at ω2. We consider representative normalized decay rates shown in the table (assuming γ = 10 −5 ω1, same as in fig.2) to explore thermal radiation features in (b,c,d). Figures (b,c) demonstrate the nonlinearity-induced changes in the thermal radiation power and statistics of intensity fluctuations g (2) (0) respectively, at frequencies ω2 (left figures) and ω1 (right figures). The inset in fig.(b) compares the maximum thermal emission in the nonlinear regime with maximum linear thermal emission. (d) Because of the nonlinear coupling, the resulting thermal emission exhibits correlations between the intensities collected at frequencies ω1 and ω2. These figures demonstrate that the strong nonlinear effects are observed when the favorable channels γ 1d , γ2e, κ (red arrows) dominate the spurious channels γ1e, γ 2d (black arrows).
thermal emitter satisfies the rate-matching condition (γ 2d = γ 2e = γ 2 /2) [16] and emits the maximum power vii given by P max 2 (0) = γ2 2 ω 2nω2,T d . This linear limit can be surpassed in the nonlinear regime where the maximum power, P max 2 (κ) = 2γ 2 ω 2nω2,T d , is obtained when the favorable channels dominate (γ 1d , κ, γ 2e γ 2d , γ 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 enhancement of frequency-selective thermal emission can be surpassed via nonlinear spectral redistribution of energy. This finding is important in the context of recent works on the far-field emission enhancements from subwavelength emitters [23][24][25] and the fundamental bounds on radiative heat transfer [28][29][30][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(c) depicts the modification in g (2) (0) as a function of nonlinear coupling κ for both resonators. Evidently, super-bunching g (2) (0) > 2 can also be observed for finite nonlinear coupling. Furthermore, figure 3(d) shows that the emitted thermal radiation collected by detectors at distinct frequencies ω j for j = [1,2] can exhibit correlations. We quantify these correlations as, C = I(ω1)I(ω2) I(ω1) I(ω2) − 1, where the correlations are zero when there is no coupling between the resonators (C = 0 when κ = 0). As demonstrated, the correlations are nonzero for finite nonlinear coupling κ 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 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][45][46][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][45][46][47]. It was shown in ref. [47] that using typical materials such as AlGaAs (χ (2) ∼ 100pm/V) in optimized microposts and grating structures, the nonlinear coupling κ = β λ1 (derivation provided in the supplement) can be realized. For a wavelength of λ 1 = 2µm., this results in the nonlinear coupling κ = 3 × 10 −8 ω. Since dissipative and radiative decay rates should be of the order of the nonlinear coupling κ to observe appreciable nonlinear effects, this will require the same resonators of quality factors of Q ∼ 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 (χ (2) ∼ 10 5 pm/V). A combination of these two approaches can potentially realize κ ∼ 10 −5 ω (as explored in fig2 and fig3), lowering the quality factor requirements to Q ∼ 10 5 to observe the predicted nonlinear effects.
We also note other 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 (χ (2) ∼ 10 5 pm/V) [50], enhanced nonlinearities in two-dimensional materials such as bilayer graphene (tunable χ (2) ∼ 10 5 pm/V) [53], high-Q bound states in the continuum [54,55] and inverse-design optimization [56]. Because of many alternatives and incessant progress, 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.

V. 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 ( ω 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 viii 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].
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][24][25]32] and bounds on radiative heat transfer [28][29][30][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) = 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.

Funding.
This We consider a doubly resonant system containing a χ (2) nonlinear medium and supporting resonances at frequencies ω1 and ω2 = 2ω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 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 [1]. 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: Here ξ jl for j = [1,2] and l = [d, e] denote the uncorrelated white noise sources corresponding to intrinsic dissipative (d) heat bath and extrinsic (e) environment respectively. They satisfy the frequency correlations, where · · · is thermodynamic ensemble average and T l is the temperature. The coefficients D j are the fluctuationdissipation (FD) relations which are to be obtained from thermodynamic considerations as we describe further below. The nonlinear coupling parameters β j for j = [1,2] can be calculated using standard perturbation theory approach for weak nonlinearities [2,3]. They depend on the overlap integral of the linear cavity fields E(ω j ) and are given below: It follows that β 1 = β * 2 = β for ω 2 = 2ω 1 (frequency matching condition). The power flows leaving the modes via nonlinear coupling are:  fig.2(a) of the main text. Evidently, the semi-classical theory leads to significant deviation from thermal equilibrium for finite nonlinear coupling and despite perfect frequency matching (ω2 = 2ω1).
By substituting a j = ω jãj for j = [1,2] in above equations so that ã * jã j = |aj | 2 ωj denotes the mean photon number and then comparing with corresponding expressions in the quantum theory (equations (7) and (8) in the main text), it follows that the nonlinear coefficient β is related to the nonlinear coupling κ through the relation By making this analogy, the nonlinear coupling κ used in the quantum theory can be computed from the overlap integrals in (4).
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 ) ii which is: with Fokker-Planck coefficients, K aj a * j = K a * j aj = D 2 j Θ(ω j , T d ) + 2γ je Θ(ω j , T e ), K aj aj = K a * j a * j = 0, cross terms K aj a k (forj = k) = 0 The above coefficients are derived based on the Ito interpretation of stochastic calculus [4,5]. In the classical regime ( ω j k B T ), we can use the Maxwell-Boltzmann probability distribution P = e −U/k B T where U = |a 1 | 2 + |a 2 | 2 is the cavity energy. Here, the nonlinear contributions to the cavity energy can be safely ignored since the nonlinearity is considered perturbatively in the derivation of coupled mode equations [1]. Using the above known probability distribution in the classical regime, we derive the FD coefficients D j = 2γ jd where γ jd is the corresponding dissipation rate. The same FD coefficients can be further justified using a microscopic classical theory of damping. The nonlinearly coupled oscillators are linearly coupled with a continuum of classical heat bath oscillators of different effective temperatures k B T 1 = Θ(ω 1 , T ) and k B T 2 = Θ(ω 2 , T ). Since the nonlinear coupling does not affect the the system-heat bath linear coupling in this microscopic theory, the resulting Langevin equations have the same FD relations as in the linear regime.
We numerically simulate the stochastic equations in (1) and (2) and compute the mean photon numbers in the steady state as n ωj = |aj (t)| 2 ωj . Figure 1 compares the semi-classical theory with quantum theory for the exact same parameters as considered in figure 2(a) of the main text. In the figure, we have used the relation given by (7) to demonstrate the dependence of mean photon numbers on the nonlinear coupling κ instead of the nonlinear coefficient β. Evidently, the steady state values obtained using semi-classical theory indicate significant deviation from equilibrium mean photon numbers of heat-bath oscillators as a function of nonlinear coupling. Since a large relative deviation from thermal equilibrium occurs despite the perturbative nature of nonlinearity (κ ≈ 10 −5 ω) and the perfect frequency matching condition ω 2 = 2ω 1 , we conclude that the semi-classical Langevin theory fails to accurately describe this particular nonlinear system. We have also undertaken a different (unsuccessful) approach of multiplicative noise where the terms D j in the coupled mode equations [(1) and (2)] are dependent on mode amplitudes a j for j = [1,2]. Their unknown functional form can be obtained using the Fokker-Planck equation by enforcing the known equilibrium probability distribution.
However, we refrain from discussing it here since the discussion of multiplicative noise further requires physical justification for their particular (complicated) mathematical form [5].
We note that both semi-classical and quantum Langevin equations agree with each other in the linear regime or for (effectively or otherwise) linear form of coupling [6][7][8]. The semi-classical Langevin equations are preferred to the master equation approach particularly when many resonant modes are involved. However, for the nonlinear form of the interaction between the resonant modes of distinct frequencies, only the quantum theory is found to be adequate. Within the quantum formalism, the solution is obtained by numerically solving the corresponding quantum master equation as discussed in the main text.