Ultra-coherent meta-emitter tailors arbitrary thermal wavefront

Domain: neuroscience
Published: 2026-01-31
Editorial outcome: minor_revision

Abstract
Conventional coherent light-field manipulation techniques inherently conflict with the spatiotemporal incoherence of thermal radiation sources. While recent advances in thermophotonics have facilitated directional thermal emission, arbitrary thermal wavefront control—a cornerstone for advanced functionalities like focusing and holography—remains an unaddressed challenge. Here, we report a generalized recipe of designing meta-emitters with lossy and lossless outer boundaries, that enables thermal emission with arbitrarily tailored wavefront. Lossy and lossless surfaces on two sides of the meta-emitter are synergistically coupled by a single-mode waveguide, transforming incoherent thermal photons to coherent surface waves for wavefront shaping functionalities. Designer surface mode of meta-emitter permits the independent optimization of photon lifetime and propagation length, thus enabling scalable spatial coherence engineering. For the proof of concept, we experimentally demonstrate near-diffraction-limited self-focusing emission, quasi-two-dimensional (quasi-2D) high-quality thermal holography without speckle noise and spatial-multiplexed holography. Coupling optimization further suggests that spatial coherence exceeding 1000λ are achievable. Our proposed meta-emitter establishes a paradigm-shifting framework to integrate stochastic thermodynamic emission with precision photonic engineering, opening avenues for information-rich thermal radiation technologies.

Introduction
As a ubiquitous physical phenomenon rooted in fundamental thermodynamics, thermal radiation inherently generates broad-spectrum electromagnetic (EM) waves governed by Planck’s law1. While its applications span from primitive illumination to a variety of modern technologies such as infrared detection2 and passive cooling3, conventional implementations predominantly manipulate photon intensity within the diffusive regime. The inherent spatiotemporal incoherence and polarization randomness of thermal photons set an insurmountable barrier for their integration into phase-critical photonic systems—a realm exclusively served by coherent laser sources. Although bulk optical configurations incorporating spectral filters and polarization modulators can partially enhance thermal light correlation, such approaches inevitably introduce prohibitive energy losses and system complexity4. This persistent challenge underscores the critical scientific question: how to engineer compact thermal emitters that directly generate wavefront-engineered coherent radiation?

Resolving this conundrum holds transformative potential for bridging the gap between thermodynamic principles and photonic information technologies, representing an enduring pursuit in the scientific community over the past decades.

Enhancing thermal coherence necessitates simultaneous phase correlation improvement across both spectral and spatial domains. In 1990s, Greffet et al. discovered long-range near-field correlation can be achieved along the emission surface with the aid of resonant surface modes, such as surface-plasmon polaritons (SPPs) or surface-phonon polaritions (SPhPs)5,6. In 2002, they pioneered the experiment to generate and observe this coherent thermal emission over a polar crystal (silicon carbide) surface7. Subsequent advances in metasurfaces spurred the development of resonant architectures—including SPhPs8–13, SPPs14–18, coupled resonant cavities19, bound states in the continuum (QBIC)4,20,21 and extended guided modes22–24—enabling the engineering of spatiotemporal coherence through tailored dispersion relations.

For these resonant absorptive surfaces, resonant modes rely on collective interactions across a large number of periodic units, rendering their resonant behavior highly sensitive to the structural variations. However, structural complexity is crucial for shaping wavefronts25. While slightly introducing local perturbations can achieve asymmetric unidirectional emission21, band splitting introduces additional lattice modes that severely disrupt the emission wavefronts and degrade the signal-to-noise ratio (SNR), as illustrated in the left panel of Fig. 1a. Further increasing the structural complexity to an aperiodic arrangement necessitates multi-layered nano-fabrications to suppress the local perturbation to the nonlocal resonances, which hinder its experimental developments4. Besides, long range of disorder leads to a higher density of split subbands, approaching a quasi-continuum band (see Supplementary Notes S1.2-1.3 for detailed analysis). This broadened band degrades the coherence of resonances, making it challenging to achieve complex thermal wavefront emission, such as thermal focusing and holography.

Separating the lossy resonance and lossless emitter could significantly suppress the local perturbations to the resonance, thereby theoretically enabling complex emission functionalities26. However, their experimental realization faces challenges, such as position-dependent absorptive-structural fabrication and limited emission power. Due to these limitations, although extensive efforts have been made to enhance the coherence of the thermal radiation, all the experimental demonstrations so far have been confined to planar wavefront emission. More complex wavefront emissions, such as thermal focusing and holography, remain largely theoretical4,13,26. Detailed emission performances and mechanism comparisons between our work and previous studies can be found in Supplementary Note S1.1 and Supplementary Table 1.

Here, we introduce a thermal meta-emitter incorporating separating but synergistic waveguide-coupled lossy and lossless boundaries, enabling arbitrary thermal wavefront manipulation capability without sacrificing SNR or too much emission power, as shown in Fig. 1a. The single connecting waveguide ensures fixed scattering phase relationship on the scattering surface, endowing meta-emitters with arbitrary phase modulation capability. Experimentally, we demonstrate self-focusing thermal emission with near-diffraction-limit performance and quasi-2D high-quality coherent thermal holography with negligible speckle noise, providing the first experimental demonstration of deterministic thermal wavefront control that approaches monochromatic-level performance metrics. Furthermore, we show that the proposed framework can be spatially encoded to realize multiplexed thermal holography. By optimizing the connecting resonance, the spatial coherence of our device can be further extended to 100λ₀ while retaining all these fascinating functionalities (asymmetric emission, focusing and holography). This work establishes a paradigm-shifting platform for coherent thermal radiation control, enabling functional wavefront shaping capabilities that bridge the gap between thermodynamic emission and precision photonic engineering.

Results
Architecture of meta-emitter
The meta-emitter incorporates near-field synergistically coupled lossy and lossless boundaries, which can be aptly conceptualized as a dual-funnel (DF) model, as schematically shown in Fig. 1b. Similar to conventional thermal metasurfaces, the bottom ‘lossy funnel’ engineers the thermally populated electromagnetic mode, and transports it to the “lossless funnel” through a coupling resonant “tunnel” while extending their lifetime. These long lifetime photons correlate the scatterers in the “lossless funnel” with designed phase-delay, thereby efficiently enabling the transformation of the temporal coherence (W t) to the spatial coherence (W x) of tailored emission waves. The connecting tunnel serves as an effective coherent source, fixing the scattering phase on the lossless funnel, similar to the point source in ref. 26. The thickness of the design can be suppressed to 1λ₀, which is similar to that of conventional thermal metasurfaces.

The corresponding physical architecture is shown in Fig. 1c, comprising two corrugated metallic surfaces with an inner bridging waveguide (WG) tunnel at the center. In this bull’s eye configuration27–30, thermally populated photons from the lossy bottom surface are efficiently transported to the lossless top surface for devising manipulation, driven by the strong Fabry-Pérot (FP) resonance occurring inside the central WG tunnel. Besides, this narrow WG tunnel effectively mitigates the dispersion-induced “rainbow effect” prevalent in periodic structures7,14,31. The cavity mode’s quality factor (Q-factor) dictates the temporal coherence of transmitted thermal photons, which is transformed into the spatial coherence through the propagation of surface modes. Local scatterers on this radiation surface can be strategically introduced as defect perturbations to contrive phase-correlated thermal emission.

Different from the previous studies usually utilizing band edge standing wave modes with a large Q-factor to enhance spatial coherence, our architecture decouples resonator and scatterers to fully leverage propagating surface wave to efficiently transform the temporal coherence into spatial coherence. The spatial coherence L of meta-emitter characterizes the distance traveled by wave packet accumulating a specific phase-delay difference, which can be linked to the temporal coherence via the relation Δϕ=LΔβ=τΔω, where Δβ and τ are respectively the wave vector difference within Δω and the photon lifetime. For a high Q resonance (Δω→0) with a specific mode lifetime, the spatial coherence can be expressed by the function of mode lifetime as L ≈ τ dω/dβ = τ * vg, where vg denotes the group velocity. Supplementary Note S2 showcases an example of the two surface modes with similar lifetime at the same operational frequency could have obviously different spatial coherences along the surface. For the band edge mode, the band variance can no longer be reflected by the group velocity; the higher-order polynomial expansion term, such as band curvature (d²/dω=dβ²), needs to be applied4,21,32 to estimate the spatial coherence.

Physically, this conclusion can be understood as the distance traveled by wave packet within its lifetime. It’s worthy to note that, temporal coherence (Q-factor of the mode) just set an upper boundary of spatial coherence (L = cτ, where c is the light speed in free space), different spatial coherence can be obtained with the same temporal coherence by modifying the slope of the operational band, which is not clarified in the previous studies.

Leveraging this temporal-to-spatial coherence transformation concept, we employ spoof SPPs in our study as energy carriers owing to their low-loss and band-engineerable characteristics across the entire thermal spectrum spanning from terahertz to mid-infrared frequencies. Critically, groove-geometries enabled radiation phase and group velocity manipulation capability makes it easily meet the versatile requirements for complex thermal wavefront shaping. Figure 2a plots the band diagrams of three representative spoof SPP modes explored in our design: a non radiative bounded mode (BM) below the light cone and two radiative modes: a weakly radiative mode (W-RM) and a strongly radiative mode (S-RM). The W-RM can be viewed as a perturbed version of the BM, where enlarging one of the grooves within a long period folds the band into the light cone. The S-RM represents the extreme case of the W-RM, obtained by removing all the interconnection grooves within the period. The band diagrams of three modes with light lines are shown in Supplementary Fig. 5. The plasma frequency follows the relation ω_spoof = πc/2h (where c is the light velocity in vacuum), allowing for broadband adaptability through freely adjustable groove depth h33,34. As shown in Fig. 2b, with realistic material loss, the surface mode lifetime τ_surf (defined as τ_surf = Q_surf/ω₀) and decay length L_D (calculated as: L_D = τ_surf * vg) of spoof SPPs can vary in two orders of magnitude, providing the flexibility in accommodating devices in different dimensional scales.

When the propagation loss is small, the spatial coherence of thermal photons is mainly influenced by the finite lifetime τ_WG of WG mode, which determines the bandwidth (Fig. 2c, inset). This WG mode exhibits a dynamic balance between radiative dissipation and intrinsic damping losses. As shown in Fig. 2c, τ_WG (Q-factor, τ_WG = Q_WG/ω₀) peaks at τ_WG,max = 25τ₀ with τ₀ = 2π/ω₀ for optimal WG width w₁, confirming the geometric configuration that minimizes loss. Our theoretical framework (Discussion and Supplementary Note S14) predicts that advanced tunneling cavity designs could suppress loss channels, potentially enhancing τ_WG,max by two orders of magnitude. The energy transfer efficiency from bottom surface mode to the radiation interface achieves resonant amplification, showing a significant enhancement at matched frequencies compared to off-resonance states (Fig. 2d). This phenomenon arises from near-field coupling between confined surface waves and the WG cavity mode, as visualized by their superimposed electric field distributions (Fig. 2e). The hybrid modes create a photonic “superhighway” that preferentially channels thermal energy while suppressing parasitic dissipation. The lifetime of the hybrid mode is governed by the reciprocal relation 1/τ = 1/τ_WG + 1/τ_surf. Detailed efficiency simulations and comparisons between our design and classical giant transmission schemes27,28,35 are provided in Supplementary Note S3.1. By applying an additional metal board to transform the surface waveguide into a photonic crystal, the bottom surface can be designed as a defect cavity, boosting the emissivity of the device to ~75%. The corresponding design strategies and discussion are also provided in Supplementary Note S3.2.

Crucially, our DF architecture enables independent modulation of these critical parameters. Figure 2g illustrates the spatial coherence of three modes (BM, W-RM, S-RM) as a function of the geometric parameter w₁ at the operational frequency ω₀. The WG mode’s finite lifetime fundamentally limit coherence, with measured maxima of 15λ₀, 14λ₀ and 9λ₀ for BM, W-RM, and S-RM, respectively. This discrepancy highlights propagation loss effects under high-Q bridging cavity conditions. For validation, cross-spectral density components W_zz(ρ, z = λ₀/10) = E_z*(O, z)E_z(ρ, z) are plotted in Fig. 2f. The coherence lengths of the W_zz are 14λ₀, 13λ₀ and 8λ₀ for BM, W-RM and S-RM, respectively, consistent with the monochromatic model predictions, laying the foundation for the further wavefront manipulation. Detailed derivations are provided in the Method section.

In the following, advanced coherent devices can be obtained through the rational design of local defect scatterers on the emission surface. Due to the narrow width of the designed scatterers (around λ₀/10), the emitted phase can be simplified to solely depend on the propagation phase of the spoof SPPs38,39. To streamline the experimental implementation and analysis, our experiments are conducted in the terahertz regime where wavelength-scale phenomena permit precise defect engineering and near-field characterization.

Self-focusing thermal emission
Our first application is to demonstrate diffraction-limited self-focusing thermal emission using a one-dimensional (1D) meta-emitter. As schematically depicted in Fig. 3a, the self-focusing device features a hybrid mode design. The bottom surface is engineered to operate in S-RM, enabling efficient harvesting of ambient thermal energy. Meanwhile, the top surface functions in W-RM and is divided into two distinct functional zones: an outer radiative region (S_r) and an inner coupling region (S_c). Grooves in the coupling region are optimized to match momentum with the WG mode, whereas those in the radiative region are engineered to operate in W-RM. Here, local scatterers are strategically introduced to fulfill the desired emission wavefront. For focusing functionality, it will satisfy a spherical spatial phase profile, i.e., φ(x) = k₀ * (√(x² + f²) - f) (1)

the condition of φ(x) + k_sspp * x = 2mπ, where m is an integer and k_sspp is the wave vector of spoof SPP optimized for radiative grooves. Details on parametric optimizations, including k_sspp, groove dimensions, and S_c functionality, are provided in Supplementary Note S4 and Supplementary Table 2. For the optimized structure, the eigenmode field analysis (Fig. 3b) reveals a well-defined focusing pattern with a focal length f = 10λ₀ (λ₀ ≈ 3mm) and numerical aperture N.A. = 0.7. It achieves a high focusing efficiency of ~95% (defined as the ratio of the integrated energy at the focal line to the energy at the WG exit, details can be found in Supplementary Note S4 and Supplementary Fig. 10) and a tightly confined focal line of 0.93λ₀ (equivalent N.A. = 0.53) realized under the incidence of nonuniform surface wave across the aperture. The focal line width can be further reduced to the ideal limit of 0.72λ₀ by engineering the radiative loss of the top surface to achieve an near-uniform scattered amplitude wavefront at the expense of dissipating surface wave energy at the edge of the device. A comprehensive time-domain analysis in Supplementary Note S5 elucidates the transient focusing process initiated by a dipole source positioned near the bottom surface. Supplementary Movie 1 vividly exhibits the dynamic evolution of the hybrid coupling of the spoof SPPs with WG mode, significantly extending the lifetime of thermal photons. In addition, it is important to note that although the FP modes in the central WG exhibit dispersion with respect to the transverse wave vector k_y, as shown in Fig. 3c (detailed information can be found in Supplementary Note S6), the leaky groove-based metalenses employed here feature the desired achromatic focusing property (Fig. 3d). Photons of different frequencies converge to form a transverse focal line at the common focal plane z = f, which boost energy utilization efficiency.

To numerically characterize the stochastic thermal emission behavior of our design, we performed 100 simulations with randomly distributed dipole excitation below the bottom surface with distance of 0.3λ₀ to obtain statistical results (detailed simulation methods given in “Method”)26. Averaged intensity profiles across different planes (Fig. 3e–g) demonstrate robust focusing performance: the obtained full-width-at-half-maximum (FWHM) (1.03λ₀) of the ensemble-averaged focal line closely approaches the diffraction limit (0.93λ₀) observed in eigenmode simulations (Fig. 3b), confirming high spatial coherence across the entire aperture. The slightly broader FWHM in the ensemble-averaged case relative to the eigenmode result (inset, Fig. 3b) arises from contributions of off-resonant frequencies near the eigenfrequency. This effect can be further mitigated by enhancing the Q-factor of the WG mode, such as tunnel geometry optimization analyzed in the Discussion part.

With the optimized design, a copper sample was fabricated using the computer numerical control (CNC) cutting technique. The transmission spectrum under normal incidence (inset, Fig. 3h) shows good agreement with simulations, though the measured transmission peak is slightly broader and blue-shifted, primarily due to fabrication tolerances. The groove fabrication tolerances are within the sculpturing precision of state-of-art CNC techniques. A weak absorptive layer, ~60 μm blackbody paint (BB paint Musou black 2.0, n≈1+0.02i @100GHz), was applied to the bottom surface for thermal photon generation, as shown in the inset of Fig. 3h. This BB layer has minor influence on the lifetime and decay length for S-RM mode in the design, due to its weak localization property. The detailed analysis on the lifetime and decay length degradation for surface modes caused by BB layer is presented in the Supplementary Note S7. The detailed fabrication tolerance analysis is presented in Supplementary Note S8.

The thermal measurement setup and sample photographs are presented in Fig. 3h. The detailed experimental information is given in the Methods section. Using this setup, intensity profiles at different planes were mapped, as shown in Fig. 3i–k. The intensity mapping validation is analyzed in Supplementary Note S9. The FWHM of the measured focal line is ~3.2mm (1.06λ₀), closely matching the statistical numerical result (1.03λ₀) plotted in Fig. 3e. Importantly, the near-diffraction-limit focusing performance indicates that all grooves contribute coherently to the focal spot, with a spatial coherence equivalent to the aperture size (11λ₀)—almost an order of magnitude longer than the state-of-the-art reported in ref. 21. Other measurements align well with simulations, including intensity distributions in xz and yz planes, as quantitatively compared in Supplementary Fig. 19c–e. The intensity local minima near the center WG is originated from the non-ideal interference of the waves scattered from the left and right sides of lens grooves, due to the depth fabrication error in the coupling region around the slit exit. The detailed analysis can be found in Supplementary Note S8.4.

The flat-topped focal line profile arises from the detector’s finite response band (75–120GHz), which attenuates high-frequency components with large transverse wave vector k_y, which is analyzed in Supplementary Note S10 by ray-tracing calculation. Detailed experimental methodologies are given in the Methods section.

Polarization and spatial multiplexed coherent thermal holography
In this section, we demonstrate 1D and quasi-2D coherent thermal holography with more complex light-field manipulation using our DF meta-emitter architecture. Leveraging the design flexibility of S-RM, we introduce polarization multiplexing along orthogonal directions. The shallow grooves (<λ₀/5) together with low filling ratio (<10%, defined as groove width w to the period p, F = w/p) of S-RM, ensuring that the cross-polarization-coupling between orthogonal channels can be neglected. By integrating these advantages with a scattering-based holographic design approach39–41, we experimentally realize quasi-2D coherent thermal holography through superposition of patterns with distinct polarization states, as exemplified by the digital “8” in Fig. 4a–d. Higher-order diffraction observed in Fig. 4b, d—attributed to fabrication tolerances—can be minimized using precision nano-fabrication techniques. Non-uniform background noise arising from detector zero-point drift is addressable via advanced methods like lock-in amplification42,43. Experimental results for additional holographic numbers (0, 4, 7) are shown in Fig. 4e–g, all achievable through the proposed design framework. To quantitatively evaluate the quality of holographic image44, we calculate the SNR = <I>/σ, in which <I> and standard deviation of intensity are the mean and standard deviation of intensity within the blue rectangle box marked in Fig. 4f. Moderate coherence of the engineered thermal light greatly suppresses the speckle noise45, leading to a high image quality (SNR = 22). Detailed design schematics and polarization multiplexing strategies are provided in Supplementary Note S11.

To further showcase the wavefront control versatility of our design, we develop and experimentally realize a spatially multiplexed holographic device, as depicted in Fig. 4h. Grooves within the central shadowed region of the top surface are optimized to simultaneously accommodate two distinct holographic wavefront profiles (Fig. 4i). By leveraging the directional propagation of spoof SPPs along the radiative surface, two unique emission states are achieved through selective excitation of corresponding WGs. Experimental switching between blocked WGs yields separate holographic patterns of ‘5’ and ‘2’ (Fig. 4j–k), with the activated excitation slit WG indicated below each image. Detailed optimization procedures and structural parameters are provided in Supplementary Note S12. Additionally, adjusting the WG width enables tuning of emitted wave coherence: off-axis focusing with multiple diffraction orders demonstrates low coherence, while spatially multiplexed designs with insufficient mode decay highlight excessive spoof SPP spatial coherence (Supplementary Note S13). This demonstration of polarization and spatial multiplexing in thermal holography underscores the design’s feasibility and adaptability for complex wavefront engineering.

Methods
Simulation method of surface mode
The dispersion properties and mode profiles of BM, W-RM and S-RM in Fig. 2a are obtained through 2D simulations using the eigenfrequency solver in COMSOL Multiphysics. The material of the metal is set as copper (σ = 5.8 × 10⁷ S/m) in the simulation. The propagation lengths along the surface for each mode in Fig. 2b are calculated from the Q factors and wave vectors of modes at their designed frequencies. The Q factor indicates how quickly the mode loses energy during one cycle of oscillation and can be directly reflected in the lifetime of the mode:

τ = Q/ω (2)

where ω₀ is the angular frequency of the resonance.

The lifetime τ characterizes the decay property in the time domain, defined as the time of the mode energy decay to 1/e of its original value. For its counterpart, the propagation length of the spoof SPPs reveals the spatial distances over which the mode energy decays to 1/e of its original value, which can be expressed as

L = vg * τ (3)

where vg is the group speed of spoof SPP along the surface, given by

vg = dω/dβ (4)

Substituting for vg and τ, we can rewrite the expression for L as:

L = (dω/dβ) * Q/ω₀ (5)

where Q factor and dω/dβ can be obtained from simulation.

The eigenfrequency, Q facor, and eigenmode profile shown in Fig. 2c are obtained by eigenfrequency study in COMSOL Multiphysics. Perfect Matched layers are set as all the boundaries.

The spatial coherence of the entire device presented in Fig. 2g is calculated by the propagation length of the hybrid mode of slit mode and surface mode on the top side. The Q factor of the hybrid mode can be obtained by:

Q_hybird = (Q_surf⁻¹ + Q_WG⁻¹)⁻¹ (6)

Then, the spatial coherence can be obtained by the Equation 5.

Calculation of the cross-spectral density component W_zz
To characterize the spatial coherence of the spoof SPPs on the top surface, we calculated the cross-spectral density component of W_zz for the electric field of spoof SPPs (E_z) above the surface at z = λ₀/10, with λ₀ = 3mm. The E_z field distribution is obtained by performing 2D frequency domain simulation in COMSOL Multiphysics software. Perfectly matched layers are set as the boundary in the simulation. To attain the E_z fields at different frequencies, an electric dipole with a frequency range from 75–120 GHz is set as the source to excite the spoof SPPs on the bottom side. Three types of top surfaces are simulated. The space-frequency dependent real parts of E_z components are drawn in Supplementary Fig. 5a–c, respectively.

The slit resonance determines the frequency-dependent part, and the same slit results in similar sharp peaks for all the surface modes. The spatial dependence of E_z reflects the energy attenuation of surface modes across the surface, which corresponds to the propagation lengths of modes shown in Fig. 2b. For S-RM, the short propagation length (~10λ₀) leads to obvious attenuation in E_z, as illustrated in Supplementary Fig. 5a. In contrast, the strong EM wave confinement capability of BM and W-RM enables them to exhibit a negligible attenuation of E_z within the range of 20λ₀, Supplementary Fig. 5c. Cross-correlation function of E_z for these three modes at a different spatial interval of ρ can be expressed as36,37:

w_zz(ρ, z = λ₀/10, ω) = E_z*(O, z, ω)E_z(ρ, z, ω) (7)

as shown in Supplementary Fig. 5d–f. The spatial attenuation behavior of E_z drawn in Supplementary Fig. 5a–c can also be reflected in the cross-correlation function at each frequency. By taking the ensemble average of w_zz, cross-spectral density component W_zz = <w_zz>_ω is obtained as described in the main text.

Simulation of the thermal self-focusing performance
We employ Ansys Lumerical FDTD Solutions for the 3D simulation of the entire design. An electric dipole with a random phase, orientation angle, and x, y coordinates is placed on the bottom surface (0.3λ₀ away from the bottom surface) as the source to capture the intensity profiles at different frequencies for each plane. This simulation is repeated 100 times, each with a different random dipole source. The average of all simulated intensity profiles is then calculated to obtain the statistical emission profiles for each plane, as shown in Fig. 3e-g. In the simulations, the emission intensities of all random dipoles are set identically. The transmitted energy obtained from each dipole simulation is shown in Supplementary Fig. 11, where the differences originate from the varying coupling efficiencies between individual random dipoles and the bottom surface mode. This excitation-efficiency variation reflects the fundamental process by which thermally populated surface waves are generated, namely that dipoles with near-normal surface polarization possess higher excitation efficiency. Such coupling-induced transmittance fluctuations do not affect the coherence enhancement performance of the designed resonant tunnel.

Measurement and characterization
To characterize the transmission property of the fabricated sample, a vector network analyzer (VNA) with thorn antenna is applied. The transmission spectrum in Fig. 3h is measured under normal incidence.

For the thermal measurement, the sample painted with a blackbody on the bottom side is heated by a hot plate. A 5mm thick metal spacer is precisely positioned between the sample and the hot plate using thermally conductive adhesive to maintain a fixed separation distance. A foil-wrapped board is applied to block background thermal noise from the hot plate. The thermal radiation pattern generated by the sample on the top side is measured by a polarization-dependent detector on a translation stage point by point. To improve the accuracy of the mapping measurement, a waveguide probe with a sharp tip is designed for the detector. The detailed information of probe design can be found in Supplementary Note S9.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript titled “Ultra-Coherent Meta-Emitter Tailors Arbitrary Thermal Wavefront” presents a dual-funnel (DF) meta-emitter architecture that achieves arbitrary thermal wavefront shaping, including self-focusing and holography, by leveraging the synergistic coupling of lossy and lossless boundaries. Although coherent thermal emissions based on nonlocal metasurfaces have been studied for a while, this is one of the few works that exploits the use of thermal emissions as a coherent source for advanced wavefront engineering applications. The concept is physically insightful, and the experimental demonstration in the Millimeter wave regime provides convincing validation.\nOverall, this study is interesting, provides valuable insights into the control and use of thermal radiations. Several technical and conceptual issues need to be clarified. Once these points are fully addressed, the manuscript can be considered for publication in Nature Commnications.\n1. In the manuscript, the following sentences are used: “These long lifetime photons correlates the scatterers in the “lossless funnel” with designed phase-delay, thereby efficiently converting the temporal coherence (Wt) to the spatial coherence (Wx) of tailored emission waves.”\n“The cavity mode’s quality factor (Q-factor) dictates the temporal coherence of transmitted thermal photons, which is converted to the spatial coherence through the propagation of surface modes. Local scatterers on this radiation surface can be strategically introduced as defect perturbations to contrive phase-correlated thermal emission.”\nIn these sentences, the expression “converting Wt into Wx” does not represent a strict conservation-type transformation. Rather, it should be regarded as an engineering picture, where long-lifetime photons are utilized to facilitate phase locking across a large-scale array. For the sake of rigor, the term “convert” is not appropriate in this context and should be replaced with a more precise wording.\n2. One of the main features of the structure is the resonant funnel, and efficiency around 10% is given in Fig. 2(d). However, important information is missing here. What is the total length of the structure in the x direction? Does this length affect the overall efficiency? Is the total input power on the lossy side taken into account? It is surpring to me of an absolute effiency around 10% since only a single funnel is used. This value is much higher than the extrodinary light transmission originally in Thomas Ebbessen’s work.\nFurthermore, the reported energy tunneling efficiency is below 10%. While this level is sufficient for proof-of-concept demonstrations, it remains far from optimal. The manuscript would benefit from a discussion on:\n- Whether this efficiency is limited by mode mismatch, metallic loss, or FP resonance detuning.\n- What strategies (e.g., mode converters, impedance matching, higher-Q cavities) could realistically improve this value toward >50% or higher.\n- Whether there is a theoretical upper limit to the tunneling efficiency in this architecture.\n3. The spatial coherence length is approximated as L ≈ τ·vg, where 1/τ = 1/τWG + 1/τsurf. For maximizing emission intensity, it is desirable for τWG and τsurf to be matched for optimal coupling, while additional scattering losses from sidewall roughness must also be considered. However, intrinsic loss imposes an upper bound on τ. A natural question is whether, under the constraint of a fixed τ, the design of vg can be exploited to further enhance L. Providing additional simulations or analysis regarding the tuning of group velocity would significantly strengthen this argument.\n4. It seems to me that there is a mess in different coherence lengths. Please make the following terms more explicit: temporal/spatial coherence, temporal/spatial coherence length. Please also compare the spatial coherence length with the one calculated by λ/Δθ, where Δθ is the angular distribution of the thermal emission at a specific wavelength.\n5. In the self-focusing section, the ensemble-averaged focal line FWHM (1.03λ0), obtained from statistical simulations of 100 random dipole excitations, is shown to be very close to the diffraction-limited value (0.93λ0) observed in the intrinsic mode simulation. However, it remains unclear how large the intensity fluctuations are across different statistical realizations. If the fluctuation is large, how can one exclude the possibility that the apparent sharp focusing is caused by partially coherent imaging?\n6. The authors implement polarization multiplexing by patterning grooves along orthogonal directions. Then, when the groove width is narrow (e.g., ~λ0/10) and the depth is large (e.g., >λ0/4), will the localized SP (or spoof-SPP) modes on the groove walls convert the Ex and Ey components into each other via corner scattering or higher-order mode coupling, resulting in p↔s cross-polarization coupling? If significant cross-coupling exists, how is the crosstalk in the reported “polarization multiplexing” estimated?\n7. Thermal focusing and holography rely on precise phase control from subwavelength grooves. However, the manuscript lacks analysis of fabrication tolerance:\n- How do random errors in groove width, depth, or position affect focusing quality or hologram fidelity?\n- What is the maximum allowable error before performance degrades significantly?\n- Are the current CNC-machined samples near this limit?\nIncluding a tolerance analysis (via simulation or perturbation theory) would greatly strengthen the practical relevance of the design.\nSomee minor suggestions:\n1. I suggest in Fig. 2(a) that only a single unit cell of structure is used to demonstrate the mode profile for the two cases of BM and W-RM. Using multiple periods might be misleading.\n2. In line 247, the authors mentioned an order of magnitude longer of the cohernece length over that in Alu’s work. I don’t this is a fair comparison because the loss associated with the spoof SP in this work is much lower.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "I have read the present manuscript on thermal emitters with great interest, and commend the authors for their creativity, but have a number of concerns that I will address below. Some are technical, some are more conceptual, with which I will start.\nMy main concern relates to the way this work is presented: supposedly outperforming, in a direct comparison as shown in Fig. 1b, thermal metasurfaces comprising a patterned emitter. However, this comparison to me is not a fair one: this work, and for example reference 26, rely on thermal point sources (in this case the funnel), the radiation of which is then manipulated via a metasurface. This is conceptually completely different from thermal metasurfaces where the entire patterned surface is emitting. In a sense, the collection of thermal emission from the backside, funneling it through a subwavelength hole, and then manipulating it, is still cascading different optical functions (generation, inducing spatial coherence, and then shaping the wavefront sequentially). This is completely different from true thermal metasurfaces, where the emitting electromagnetic mode contains fluctuating sources that populate it, resulting in a single step emission process of designer functionality. I therefore strongly discourage the authors to make such strong claims and by making a direct comparison with “standard” thermal metasurfaces, particularly because these results are interesting, but claiming this much discredits the work.\nFurthermore, very strong claims are made about creating holograms, but in this case the holograms appear to simply be stitched-together lenses as in the first example, given that the holograms comprise only straight lines and the sample appears to be a number of combined lenses even in the photo. Again, this is an interesting demonstration, but this is not a hologram in the traditional sense, where a continuous phase profile yields a structure a certain distance away from the emitter that can comprise more than only straight lines (two cylindrical lenses next to each other also are not typically considered a hologram).\nSome technical concerns:\n- One big practical difference is that thermal metasurfaces can reach unity emissivity for specific polarization states, but in this case the subwavelength funnel fundamentally limits the amount of power that can be emitted (a single mode waveguide is also bound by a blackbody limit, it cannot transmit more power than that). Thermal metasurfaces are already relatively weak emitters even at unity emissivity (compared to coherent sources), so I think this is important to address, especially if you have only one funnel per many wavelengths: how does this affect the emitted power?\n- How is the blackbody modeled and taken into account in determining the efficiency and size of the collection area? As far as I can tell this is not mentioned, which makes me think the black body is not incorporated directly – however, if so, that would result in a misleading collection size.\n- How does this extend to the near- or mid-infrared? Thermal radiation at 100 GHz is weak, and thermal sources are sought after in the mid-IR. This is an interesting demonstration but has little practical implication: at high frequencies, one would require a very lossless material combined with a single point thermal emitter, similar to the aforementioned reference.\n- What is the spot size for the NA of the lens? It apparently is near-diffraction limited, but the value that one would expect in the ideal case is not mentioned as far as I can tell.\nTo conclude, this is an interesting demonstration of an idea very similar to ref. 26, but it is pitched in a way that is far too strong: it is a compactified set of cascaded elements (thermal source, pinhole, wavefront shaper) with very limited demonstrations of wavefront shaping (as far as I can tell only focusing has been shown, including in a demonstration of a “hologram” formed by multiple lenses), which is not what a thermal metasurface does, and I therefore don’t think it should be pitched as such or compared with them.", "reviewer_4": "The authors discuss and experimentally demonstrate a bilayer device which performs sophisticated tailoring of the wavefront of thermal emission in the millimeter-wave range. The working principle relies on separately engineering two cascaded layers, coupled via a narrow channel, to achieve more flexibility in simultaneously collecting thermal emission efficiently and radiate it with a desired wavefront.\nThe approach proposed by the authors is interesting, and, to the best of my knowledge, is conceptually novel. I can recommend publication after the authors have addressed the following comments:\n1. As a general comment, there are several information, concepts and definitions which are relegated to the supplementary material, but which I think should be mentioned in the main text in order to provide a clearer description (note that the SM is 37 pages long). I provide some example in the next comments.\n2. Can the authors provide a clear definition for the parameters tau_surf and L_D introduced at lines 141-142?\n3. Can the authors clarify what is the difference between the designs BM and W-RM. Based on the sketches in the inset in Fig. 2a, they lookthe same.\n4. Can the authors add a light-cone line to Fig. 2a, to better show which modes are dark?\n5. Line 152, “The energy transfer efficiency from source to radiation interface achieves resonant amplification, showing a significant enhancement at matched frequencies compared to off-resonance states (Fig. 2d).”.\nIt is not clear how such “transfer efficiency” is defined, and what is the “source to radiation interface”. Can the authors provide a more detailed description of how the plots in Fig.2d are calculated?\n6. At line 148, there is a typo in the definition of tau_WG=Q*omega_0. I assume the correct version is Q/omega0\n7. I am not sure I fully understand the calculations in Figs. 2f-g. Are these calculations done for a bilayer device, or for a single layer? If my understanding is correct, here the authors calculate the eigenmode of the structure (either single layer or bilayer, and for a fixed design type), and then evaluate the function W_zz for a fixed value of z=\\lambda0/10. Can the authors confirm if it is a bilayer or single-layer simulation, and if it’s an eigenmode calculation?\nIntuitively, I am not sure I understand why the spatial coherence should show such a strong dependence on the parameter W1, since W1 affects only the behavior near the center. Can the authors elaborate more on this?\n8. For the calculations shown in Fig. 3b, the authors need to add more detail the main text, such as: materials, material parameters, geometrical parameters, and value of lambda0\n9. Line 207, “It achieves a high energy focusing efficiency of ~95%”. Please clarify how focusing efficiency is defined here\n10. For the simulations in Fig. 3e, the authors say that “we performed 100 simulations with randomly distributed dipole excitation to obtain statistical results”. Can the authors clarify in which region of space the dipoles are placed? Are the dipoles only placed in the region below the bull’s eye structure? If yes, is this a realistic assumption? In a realistic scenario, any part of the device will emit thermal radiation.\n11. I notice some difference between experiments (Fig. 3i) and simulations (Fig. 3e). In particular, in the experimental data there is an intensity local minima near the center of the structure (close to the WG element), while the simulations predict a high intensity. Can the authors comment on this?\n12. When describing the fabricated structure in the caption of Fig. 3h, the authors mention that the bottom side is “coated with a blackbody layer”. The need for this layer is not mentioned elsewhere in the paper, and it’s not clear which role it plays. Is the presence of this layer captured in the simulations? Please include a more detailed description of it"}, "author_response": "Thanks for your positive evaluation and valuable comments. The point-by-point responses are listed below.\n1. In the manuscript, the following sentences are used: “These long lifetime photons correlates the scatterers in the “lossless funnel” with designed phase-delay, thereby efficiently converting the temporal coherence (Wt) to the spatial coherence (Wx) of tailored emission waves.”\n“The cavity mode’s quality factor (Q-factor) dictates the temporal coherence of transmitted thermal photons, which is converted to the spatial coherence through the propagation of surface modes. Local scatterers on this radiation surface can be strategically introduced as defect perturbations to contrive phase-correlated thermal emission.”\nIn these sentences, the expression “converting Wt into Wx” does not represent a strict conservation-type transformation. Rather, it should be regarded as an engineering picture, where long-lifetime photons are utilized to facilitate phase locking across a large-scale array. For the sake of rigor, the term “convert” is not appropriate in this context and should be replaced with a more precise wording.\nResponse: Thank you for your insightful comments and rigorous corrections. We agree that the temporal and spatial coherence are not related through a conservation-type transformation, and the terms ‘convert’ or ‘conversion’ are not sufficiently rigorous in this context. To avoid potential misunderstandings, we have replaced ‘convert’ with ‘transform’ to more accurately describe the relationship between spatial and temporal coherence.\nCorresponding Modification:\n1. In the Abstract, lines 27-28, we changed the word 'converting’ into ‘transforming’.\n2. In page 3, paragraph 1, lines 120-123, we modified ‘converting’ to ‘transformation’:\n“These long-lifetime photons correlate the scatterers in the “lossless funnel” with\ndesigned phase delay, thereby efficiently enabling the transformation of temporal\ncoherence (W) to the spatial coherence (W ) of tailored emission waves.”\nt x\n3. In page 4, paragraph 1, lines 133-138, we modified ‘converted’ to ‘transformed’:\n“The cavity mode’s quality factor (Q-factor) dictates the temporal coherence of\ntransmitted thermal photons, which is transformed into spatial coherence through the\npropagation of surface modes. Local scatterers on this radiation surface can be\n1\nstrategically introduced as defect perturbations to contrive phase-correlated thermal\nemission.”\n2. One of the main features of the structure is the resonant funnel, and efficiency around 10% is given\nin Fig. 2(d). However, important information is missing here. What is the total length of the structure\nin the x direction? Does this length affect the overall efficiency? Is the total input power on the lossy\nside taken into account? It is surprising to me of an absolute efficiency around 10% since only a single\nfunnel is used. This value is much higher than the extraordinary light transmission originally in\nThomas Ebbessen’s work.\nFurthermore, the reported energy tunneling efficiency is below 10%. While this level is sufficient for\nproof-of-concept demonstrations, it remains far from optimal. The manuscript would benefit from a\ndiscussion on:\n- Whether this efficiency is limited by mode mismatch, metallic loss, or FP resonance detuning.\n- What strategies (e.g., mode converters, impedance matching, higher-Q cavities) could realistically\nimprove this value toward >50% or higher.\n- Whether there is a theoretical upper limit to the tunneling efficiency in this architecture.\nResponse: Thank you for your question. The definition and calculation of efficiency were\npresented in the Extended Method Section of the original manuscript. Although our design shares\nstructural similarities with those employed in extraordinary optical transmission (EOT) systems, it\noperates in the near field, in contrast to the far-field nature of EOT. Consequently, the efficiency\nmetrics differ fundamentally due to the distinct excitation conditions involved. In EOT systems, a\ncoherent plane wave source is applied. The transmission efficiency can be directly defined as the\nratio of transmitted power to the incident power (𝜂 (cid:3404) 𝑃 /𝑃 ), measurable via a Vector\n(cid:2930)(cid:2928)(cid:2911)(cid:2924)(cid:2929) (cid:2919)(cid:2924)(cid:2913)(cid:2919)\nNetwork Analyzer (VNA) for the microwave band or a spectrometer in the optical band. In our\ncase, however, the excitation originates from stochastic thermal fluctuations. The total transmitted\nenergy strongly depends on the dipole position and polarization, making direct quantification\ndifficult. Therefore, we used the surface mode as the reference source and defined the coupling\nefficiency as the ratio of the energy carried by the top surface mode after tunneling through the slit\nto that of the bottom surface mode. This measures the mode-to-mode transfer efficiency through\nthe three coupled components: bottom surface waveguide, slit waveguide, and top radiation\nsurface.\nWe performed simulations in COMSOL Multiphysics as illustrated in Fig. S6a (reproduced as\nFig. R1 shown below). The size of the simulated structure is ~46λ long. A dipole source is placed\n0\nat one edge of the bottom surface (~22λ from the slit) to excite the bottom surface mode. The\n0\nincident energy (P ) from the bottom surface mode is calculated by integrating energy ~20λ away\nin 0\nfrom the excitation dipole, where the slit is removed. For the full structure, three top surface\nconfigurations (Bottom: BM; top: BM, W-RM, S-RM) are simulated with a slit WG (width:\nλ /150). Two energy integration lines (P = P + P ) are placed 2λ from the slit to avoid\n0 out o1 o2 0\ncollecting the scattered waves from the slit. The energy tunneling efficiency is calculated as η =\nP /P . The electric field distributions at the eigenfrequency of slit waveguide for these cases are\nout in\nshown in Fig. S6b-e (Fig. R1b-e). So, in this case, the total length of the structure will not affect\nthis transfer efficiency. Since only the surface mode energy is considered, radiation emitted by\ndipoles that do not couple to the surface wave is intentionally excluded. For the lossy grating (S-\nefficiency in Fig. 2d.\nFig. R1. Calculation of energy tunneling efficiency. a. Sketch of calculation method. A dipole\nsource is placed at the edge of the bottom structure to excite the bottom surface mode. Two\nintegration lines are set on the top surface near the center slit WG. The distance between the\nintegration line and the slit is 2λ . The orange dashed line is the integration line for the incident\n0\nenergy of the bottom surface mode. The slit is not constructed for the incident energy simulation.\nThe length of the integration line is λ . The energy tunneling efficiency is calculated as: η = P /\n0 out\nP . b. Electric field of the bottom surface mode at the eigenfrequency of slit WG. c-e. Electric\nin\nfield at the eigenfrequency of the slit WG for the whole structure with the structure of BM\n(bounded mode), W-RM, and S-RM (weak/strong radiation mode) on the top surface, respectively.\nThe structure of BM is set as the bottom surface for all these simulations. Black lines on the\nbottom in b-d are the energy integration lines. Scale bar: λ .\nframework relies on surface waves excited by the grooves around the slit to enhance transmission.\nBecause periodic boundary conditions are not suitable for the whole structure (a single slit with\nperiodic grooves), the transmission of the entire device strongly depends on the simulation size.\nFor this reason, most previous works have only reported relative enhancement compared to an\nisolated slit, rather than absolute transmission.\nNevertheless, approximate values can be inferred from reported results. For example, G. W. ’t\nHooft, et. al. reported their optimization could squeeze energy through the slit from around 20% of\nthe width of the entire pattern (λ = 800 nm, groove number = 11)1, implying an overall\ntransmission efficiency of ~20% for the whole device. In our design, the simulated plane-wave\ntransmission is ~14%, and the experimental value is ~10% (inset of Fig.3h), for a lateral size of\n~20λ (6 cm).\n0\noptimized the distance between the first groove and slit for both sides and designed the coupling\nregion on the top side for mode S-RM, as presented in Fig. R1d. For the reason of the low\nefficiency and how to improve it, that’s an open question, and we are happy to briefly discuss it\nwith temporal coupled mode theory (TCMT)2,3.\nCase 1: The traditional giant transmission: plane wave input and plane wave output.\nIn this case, the system can be viewed as a two-port system, and the center slit functions\nas a resonant cavity, while the grooves on both surfaces function as a mode adapter between the\nplane wave and the slit mode, as shown in Fig. R2a.\nThe temporal dynamic evolution of the mode is described by the equation:\n(cid:3031)(cid:3028)(cid:4666)(cid:3047)(cid:4667) (cid:3404)(cid:4666)𝑖𝜔 (cid:3398)𝛾 (cid:3398)𝛾 (cid:3398)𝛾 (cid:4667)𝑎(cid:4666)𝑡(cid:4667)(cid:3397)K(cid:3021)|𝑆 ⟩, (ER.1)\n(cid:2868) (cid:2868) (cid:2869) (cid:2870) (cid:2878)\n(cid:3031)(cid:3047)\n|𝑠 ⟩(cid:3404)𝐶|𝑠 ⟩(cid:3397)D𝑎(cid:4666)𝑡(cid:4667), (ER.2)\n(cid:2879) (cid:2878)\nwhere a(t) is the complex modal amplitude as a function of time, ω is the mode eigenfrequency,\n0\nγ is the nonradiative decay rate, γ , γ are the radiative decay rates to the two ports, K = [κ , κ ]T\n0 1 2 1 2\nand D= [d , d ]T are the coupling matrices between resonance and the two ports, |s > = [s , s ]T\n1 2 + 1+ 2+\nand |s> = [s , s ]T are the input and output waves.\n- 1- 2-\nBased on energy conservation, reciprocity, and time reversal invariance, we can derive\nthe relationships between these coupling coefficients:\n𝑘 (cid:3404)𝑑 (cid:3404)(cid:3493)2𝛾 , 𝑘 (cid:3404)𝑑 (cid:3404)(cid:3493)2𝛾 , (ER.3)\n(cid:2869) (cid:2869) (cid:2869) (cid:2870) (cid:2870) (cid:2870)\n(cid:3398)1 0\n𝐶 (cid:3404) (cid:4674) (cid:4675). (ER.4)\n0 (cid:3398)1\nWith time-harmonic input [exp(iωt)], the solution of ER.1 can be written as:\n𝑎(cid:4666)𝑡(cid:4667)(cid:3404)\n(cid:2895)(cid:3269)|(cid:3020)(cid:3126) ⟩\n, (ER.5)\n(cid:3036)(cid:4666)(cid:3104)(cid:2879)(cid:3104)(cid:3116) (cid:4667)(cid:2878)(cid:3082)(cid:3116)(cid:2878)(cid:3082)(cid:3117)(cid:2878)(cid:3082)(cid:3118)\nWhen we assume the input as |𝑠 ⟩(cid:3404)(cid:4670)1,0(cid:4671)(cid:3021), the transmission coefficient can be expressed as:\n(cid:2878)\nt (cid:3404)\n(cid:3046)(cid:3118)(cid:3127)\n(cid:3404)𝜅 𝑎(cid:4666)𝑡(cid:4667)(cid:3404)\n(cid:2870)√(cid:3082)(cid:3117)(cid:3082)(cid:3118)\n, (ER.6)\n(cid:3046)(cid:3117)(cid:3126) (cid:2870) (cid:3036)(cid:4666)(cid:3104)(cid:2879)(cid:3104)(cid:3116) (cid:4667)(cid:2878)(cid:3082)(cid:3116)(cid:2878)(cid:3082)(cid:3117)(cid:2878)(cid:3082)(cid:3118)\nThe transmission efficiency reaches a maximum under the following conditions:\na. Resonance condition: 𝜔 (cid:3404)𝜔 ;\n(cid:2868)\nb. No nonradiative loss: 𝛾 (cid:3404)0;\n(cid:2868)\nc. Critical coupling condition: 𝛾 (cid:3404)𝛾 .\n(cid:2869) (cid:2870)\nUnder these conditions, the transmission efficiency could reach |t|2 = 100%.\nThe reason for the relatively low transmission efficiency in previous work at resonance frequency\n(~20% in paper1) came from the mode mismatch and metallic loss. The metallic loss is obvious in\nthe optical band. Besides, the mode mismatch is the main reason for the limited efficiency. For the\nperiodic groove grating with a single slit (bottom surface in Fig. R2b), the scattering field of the\neigenmode is not an exact plane wave. Additionally, the eigenmode profile of the whole structure\ncontains the surface wave at the edge of the structure (encircled in Fig. R2b), which makes the\nexcitation efficiency degrade with only plane wave excitation.\nTo verify the intrinsic limit, we excited the system with the reciprocal eigenmode profile at the\nfocal plane and set the metal as a perfect electric conductor. Under this ideal matched excitation,\nthe transmission reached 73.8%, demonstrating that >50% efficiency is achievable in principle\n(Fig. R2c).To further enhance the transmission to reach the theoretical maximum, coupling\n4\nbetween the surface wave on both sides of the slit and the coupling between the surface wave and\nplane wave should also be further optimized.\nFig. R2. Transmission efficiency analysis. a. Sketch of the design coupled with two free-space\nports. b. The |E| field of the eigenmode of the design. The dotted circle and the orange arrow\nindicate the surface wave at the edge of the structure. c. The |E| field profile of the design excited\nwith the reciprocal field of the eigenmode at the focal plane. Input and output ports are marked in\nthe figure.\nCase 2: Surface wave input and output (energy transmission efficiency in our paper).\nDue to the incoherent nature of thermal excitations, the global eigenmode of the bottom\nsurface cannot be coherently excited by incoherent thermal dipoles at different positions. In this\ncase, the bottom surface is treated as a waveguide for the thermally excited surface waves.\nSimilar to the energy calculation in Fig. R1, we further monitored the energy transporting\nefficiency at different ports (Port 1-6, P1-P6) as indicated in Fig. R3a. We found that most of the\nenergy (75%) is reflected in P1 and only ~13% energy is transported to the top side. This low\nefficiency is fundamentally caused by the large momentum mismatch between the bottom surface\nwave (primarily k component) and the slit waveguide mode (k component). Since the slit and\nx z\nsurface waveguide are orthogonal, this mismatch is difficult to eliminate, and groove–slit spacing\nprovides limited tuning freedom, which was also applied in the optimization of the structure.\nIf device footprint constraints are relaxed, this energy efficiency can be greatly boosted\nby designing the bottom surface as a cavity. When we apply an additional constraint on the\nsurface mode, such as adding a metal board above the grooves with a specific gap, the surface\nwave mode transforms into a different guided mode with a specific bandgap, as shown in Fig. R3b.\nThe bottom surface cavity can be formed by introducing a defect in the photonic crystal. In this\ncase, Bragg reflection of the photonic crystal could further confine the thermal energy distribution\naround the slit and enhance the energy efficiency. We applied a similar simulation method shown\nin Fig. R2c, the system can reach 75% absorptivity (reciprocally emissivity) at eigenfrequency by\napplying suitable absorption in the bottom photonic crystal defect cavity. However, the photon\ncollection area is limited by the effective mode area in this case.\nat each port. P -P are the incident port, direct transmit port, backscattering port, transmit port for\n1 6\nthe left side, transmit port for the right side, and the free-space scattering port, respectively. The\narrows indicate the Poynting vector. The energy received by each port is marked in the figure. b.\nBand-diagram of the photonic crystal unit cell. Inset shows the Ey field of the eigenmode profile. t,\nh, p, and g are 1.4 mm, 0.5 mm, 2.8 mm, and 0.2 mm, respectively. c. The |E| field profile of the\ndesign excited with the reciprocal field of the eigenmode at the focal plane. Input port and\nabsorptivity A (reciprocally, emissivity ε) are marked in the figure.\nCorresponding Modification:\nWe merged the energy tunneling efficiency calculation in the Method section and the above\nefficiency optimization discussion into a new supplementary section (S3) titled ‘Energy tunneling\nefficiency calculation and optimization strategy discussion’.\nIn page 4, paragraph 2, lines 189-191, we clarified the efficiency description as follows: ‘The\nenergy transfer efficiency from the bottom surface mode to the radiation interface\nachieves resonant amplification, showing a significant enhancement at matched\nfrequencies compared to off-resonance states.’\nWe also added a brief description of the corresponding supplementary part at the end of this\nparagraph (page 4, paragraph 2, lines 197-203) as follows: ‘Detailed efficiency simulations and\ncomparisons between our design and classical extraordinary transmission\nschemes27,28,35 are provided in Supplementary Text S3.1. By applying an additional metal\nboard to transform the surface waveguide into a photonic crystal, the bottom surface can\nbe designed as a defect cavity, boosting emissivity to ~75%. The corresponding design\nstrategies and discussion are also provided in Supplementary Text S3.2.\n3. The spatial coherence length is approximated as L ≈ τꞏvg, where 1/τ = 1/τWG + 1/τsurf. For\nmaximizing emission intensity, it is desirable for τWG and τsurf to be matched for optimal coupling,\nwhile additional scattering losses from sidewall roughness must also be considered. However,\nintrinsic loss imposes an upper bound on τ. A natural question is whether, under the constraint of a fixed τ,\nthe design of vg can be exploited to further enhance L. Providing additional simulations or analysis\nregarding the tuning of group velocity would significantly strengthen this argument.\nResponse: Thanks for your insightful question. It’s a critical question in thermal wavefront\nshaping and one of the advancements of our scheme over the nonlocal thermal metasurfaces in the\nprevious studies.\nTemporal and spatial coherence reflect the correlation behavior of light vibrations in the temporal\nand spatial domains, respectively. These two measures are correlated with each other when the\nlight is emitted by a specific resonance. Importantly, different spatial coherences can be obtained\neven when the temporal coherence is fixed, simply by tailoring the slop of the operational band.\nOne representative extreme case is the comparison between our work and previously reported\nnonlocal thermal metasurfaces operating at the band edge. In Ref. 4, nonlocal thermal\nmetasurfaces with a Q factor of 128 achieved a spatial coherence of only 1.6λ . In this case, band\n0\nvariation within a narrow spectral width can not longer be captured by the band slop (group\nvelocity) any more, higher order of polynomial terms, such as band curvature (𝑏 ∝ 𝑑(cid:2870)ω/𝑑𝛽(cid:2870)),\nneed to be applied. In contrast, our design, despite having a more moderate Q factor of ~30 (the\n6\ninset of Fig. 3h in the Main Text), achieves a spatial coherence of ~11λ which is almost one order\n0,\nof magnitude higher than their work.\nTo further illustrate the role of group velocity, we designed two output surface modes with\ncomparable temporal coherence but different group velocities at the operational frequency, as\nshown in Fig. R4a. The cross-spectral density components W (ρ) of the two modes as a function\nzz\nof ρ at a fixed height z = λ /10 above the surface are plotted in Fig. R4b. Although Mode 2 has a\n0 0\nslightly larger temporal coherence than Mode 1, the W (ρ) attenuates faster than Mode 1. The\nzz\n1/√𝑒 correlation lengths of the W\nzz\nare around 42λ₀, 26λ₀ for Mode 1and Mode 2, respectively,\nwhich is close to the τꞏv expectation (40λ₀, 28λ₀ for Mode 1and Mode 2, respectively). The\ng\nslightly variation between the full-wave simulation and theoretical analysis comes from the near-\nfield fluctuation of the surface wave.\nFig. R4 Spatial coherence of surface modes with similar lifetimes but different band slopes. a.\nDispersion properties of two surface modes. The |E|-profiles of eigenmodes are shown in the right\npanel. The corresponding group velocities and lifetimes are marked in the plot. b. Cross-spectral\ndensity component W (ρ) of the surface wave as a function of ρ at a fixed height z = λ /10 above\nzz 0 0\nthe surface.\nCorresponding Modification:\nIn page 3, paragraph 1, lines 139-160, we modified the analysis of group velocity in determining\nspatial coherence as follows: “\nDifferent from the previous studies usually utilizing band edge standing wave modes\nwith a large Q-factor to enhance spatial coherence, our architecture decouples resonator\nand scatterers to fully leverage propagating surface wave to efficiently transform the\ntemporal coherence into spatial coherence. The spatial coherence L of meta-emitter\ncharacterizes the distance traveled by the wave packet accumulating a specific phase-\ndelay difference, which can be linked to the temporal coherence via the relation, where\nand are respectively the wavevector difference within and the photon lifetime. For a high\nQ resonance (→ 0) with a specific mode lifetime, the spatial coherence can be expressed\nby the function of mode lifetime as, where vg denotes the group velocity. Supplementary\nText S2 showcases an example of the two surface modes with similar lifetime at the\nsame operational frequency could have obviously different spatial coherences along the\nsurface. For the bandedge mode, the band variance can no longer be reflected by the\ngroup velocity, the higher-order polynomial expansion term, such as band curvature (𝑏 ∝\nd(cid:2870)ω/𝑑𝛽(cid:2870)),\nneeds to be applied1,6,7 to estimate the spatial coherence.\nPhysically, this conclusion can be understood as the distance traveled by wave packet\nwithin its lifetime. It’s worth to note that, temporal coherence (Q-factor of the mode) just\nset an upper boundary of spatial coherence (L = cτ, where c is the light speed in free\nc\n7\nspace), different spatial coherence can be obtained with the same temporal coherence\nby modifying the slope of the operational band, which is not clarified in the previous\nstudies.\nBesides, we added the discussed example of the spatial coherence difference of surface modes\nwith similar lifetimes but different band slopes in Supplementary Text S2.\n4. It seems to me that there is a mess in different coherence lengths. Please make the following terms\nmore explicit: temporal/spatial coherence, temporal/spatial coherence length. Please also compare the\nspatial coherence length with the one calculated by λ/Δθ, where Δθ is the angular distribution of the\nthermal emission at a specific wavelength.\nResponse: Thanks for your comments. In the original manuscript, the terms ‘spatial coherence’\nand ‘spatial coherence length’ were referred to as the same quantity. To avoid ambiguity, we\nunified all terminology to ‘spatial coherence’.\nAs discussed in our previous response, spatial coherence quantifies the correlation of optical field\nvibrations over a certain spatial separation. It can be reflected by the angular divergence property\n(λ/Δθ) for plane wave emission or other formats for more sophisticated emission wavefronts, such\nas the diffraction-limited Airy disk for parabolic focusing wavefronts.\nIn our study, the full width of half maximum (FWHM) of the focal spot serves as a practical\nindicator of the spatial coherence (effective lens aperture). To fully utilize the surface wave energy,\nwe optimized the surface wave decay to reach 0 at the edge of the grating region, which results in\nan effective aperture smaller than the physical device size. In the simulation in Fig. 3b, FWHM is\n~ 0.93λ , which corresponds to a spatial coherence of ~13λ . The experimental FWHM of the\n0 0\nfocal line is slightly degraded to ~1.06λ , the corresponding spatial coherence is 𝐿 (cid:3406)11λ . By\n0 (cid:3030) 0\nsimply engineering the radiative loss of the top radiation surface to make the energy of the surface\nwave decay to 1/e at the boundary of the device, the spatial coherence can be extended to the size\nof the device, achieving a narrower focal line.\nCorresponding Modification:\nWe unified all the terms of ‘coherence length’ to spatial coherence.\nIn page 6, paragraph 2, lines 247-250. we modified the description of the spatial coherence of self-\nfocusing device as follows: ‘The focal-line width can be further reduced to the ideal limit of\n0.72λ by engineering the radiative loss of the top surface to achieve a nearly uniform\n0\nscattered-amplitude wavefront, at the expense of dissipating surface-wave energy at the\nedge of the device.’\n5. In the self-focusing section, the ensemble-averaged focal line FWHM (1.03λ0), obtained from\nstatistical simulations of 100 random dipole excitations, is shown to be very close to the diffraction-\nlimited value (0.93λ0) observed in the intrinsic mode simulation. However, it remains unclear how\nlarge the intensity fluctuations are across different statistical realizations. If the fluctuation is large,\nhow can one exclude the possibility that the apparent sharp focusing is caused by partially coherent\nimaging?\ncentrally located resonant tunnel. All dipole-excited surface waves on the bottom surface are\ncoupled into this central tunnel, where their lifetimes are elongated by the resonance. In our\nsimulations, the emission intensities of all randomly placed dipoles are kept identical. The\ntransmitted energy for each dipole simulation is shown in Fig. R5, and the observed variations\narise solely from the different coupling efficiencies between each random dipole and the bottom\nsurface mode.\nThis variation in excitation efficiency reflects a fundamental physical process governing how\nthermally populated surface waves are generated: dipoles with a near-normal surface polarization\nexhibit higher excitation efficiency of the surface wave. Importantly, such coupling-induced\ntransmittance fluctuations do not affect the coherent-enhancement functionality of the designed\nresonant tunnel. The experimental measurements at each observation plane further validate the\nreliability of our simulation results. Therefore, the coherent-enhancement mechanism and the\nthermal-focusing performance predicted by the random-dipole excitation simulations remain\naccurate and physically consistent, and the observed coupling-efficiency–induced fluctuations are\nfully reasonable.\nFig. R5 Transmittance fluctuation in 100 random dipole excitations.\nCorresponding Modification:\nWe added the discussion of transmittance fluctuations observed in the simulations to the Method\nsection (page 12, lines 480-488), as follows: ‘In the simulations, the emission intensities of all\nrandom dipoles are set identically. The transmitted energy obtained from each dipole\nsimulation is shown in Extended Data Fig. E2, where the differences originate from the\nvarying coupling efficiencies between individual random dipoles and the bottom surface\nmode. This excitation-efficiency variation reflects the fundamental process by which\nthermally populated surface waves are generated, namely that dipoles with near-normal\nsurface polarization possess higher excitation efficiency. Such coupling-induced\ntransmittance fluctuations do not affect the coherence-enhancement performance of the\ndesigned resonant tunnel.\nFig. E2. Transmittance fluctuations in 100 random dipole excitations.’\n9\n6. The authors implement polarization multiplexing by patterning grooves along orthogonal directions.\nThen, when the groove width is narrow (e.g., ~λ0/10) and the depth is large (e.g., >λ0/4), will the\nlocalized SP (or spoof-SPP) modes on the groove walls convert the Ex and Ey components into each\nother via corner scattering or higher-order mode coupling, resulting in p↔s cross-polarization\ncoupling? If significant cross-coupling exists, how is the crosstalk in the reported “polarization\nmultiplexing” estimated?\nResponse: Thanks for your question. It’s a very important issue when we designing the\npolarization multiplexing devices. In our design, the grooves are relatively shallow (<λ /5, the\n0\ndetailed parameters are listed in the Supplementary Table S10), and the localized surface-plasmon\nresonance is not strong enough to confine substantial energy inside the grooves. As a result, most\nof the surface-wave energy propagates through the scatterers rather than being trapped in the\ngrooves.\nOn the other hand, our polarization-multiplexing design employs the S-RM mode, which does not\nrequire inter-coupling gratings between scatterers. The small filling ratio of the grooves (~10%)\nalso minimizes perturbation, further suppressing potential crosstalk between orthogonal\npolarization channels. In our experiments, the designed patterns are clearly observed, indicating\nthat the cross-polarization coupling in our design is negligible.\nCorresponding Modification:\nIn page 8, lines 315-319, we added the following analysis on cross-polarization-coupling:\n‘Leveraging the design flexibility of S-RM, we introduce polarization multiplexing along\northogonal directions. The shallow grooves (<λ /5) together with the low filling ratio\n0\n(<10%, defined as groove width w to the period p, F = w/p) of S-RM, ensuring that cross-\npolarization coupling between orthogonal channels can be neglected.’\n7. Thermal focusing and holography rely on precise phase control from subwavelength grooves.\nHowever, the manuscript lacks analysis of fabrication tolerance:\n- How do random errors in groove width, depth, or position affect focusing quality or hologram fidelity?\n- What is the maximum allowable error before performance degrades significantly?\n- Are the current CNC-machined samples near this limit?\nIncluding a tolerance analysis (via simulation or perturbation theory) would greatly strengthen the\npractical relevance of the design.\nResponse: Thanks for your questions. The tolerance analysis is listed below:\n1. Impact of geometric fabrication errors.\nDeviations in groove width, depth, and position lead to variations in the radiative phase of surface-\nwave emission, which degrade the overall wavefront modulation performance. Due to the narrow\nwidth of the designed scatterers (around λ /10), the emitted phase can be simplified to solely\n0\ndepend on the propagation phase (𝜑 (cid:3404) 𝑘 𝑥) of the spoof SPPs. Errors in groove width and\n(cid:2929)(cid:2929)(cid:2926)(cid:2926)\ndepth mainly cause the wavevector (𝑘 ) deviation, and the groove position (x) errors directly\n(cid:2929)(cid:2929)(cid:2926)(cid:2926)\ncause the propagation phase deviation. We conducted the simulation of wavevector deviation (δk)\non groove width and depth deviation (δw and δh) in one unit cell, as presented in Fig. R6a-b.\nAround 1% fabrication error in width and depth causes 1.5% and 6% of wavevector deviation (δk).\nPositional deviation (δx) approximately proportional to the scattered phase deviation (δ𝜑).\nwidth and depth variation. c. Scattered phase deviation as a function of groove position deviation.\n2. Tolerance criterion based on optical wavefront aberration.\nIn traditional optical systems, wavefront deviation (wavefront aberration) less than 1/4λ are\n0\ngenerally considered acceptable. Borrowing this criterion, the allowable fabrication tolerances for\nour wavefront-modulation device are approximately: groove width, depth, and position tolerance\nare about ±17%, ±4% and ±25%.\n3. Compatibility with CNC fabrication precision.\nState-of-the-art CNC machining can achieve groove-width, groove-depth, and position accuracies\nof approximately ±0.02 mm, ±0.01 mm, and ±0.004 mm, respectively, corresponding to\nfabrication errors of ~±10%, ±9%, and <±0.5%. These values fall within the individual tolerance\nlimits noted above.\nHowever, the accumulated phase error can become significant if all parameters deviate\nsimultaneously near their upper tolerance limits. In such a worst-case scenario, the combined\nphase error can approach ~20%, close to the allowable tolerance boundary. This accumulated\nphase deviation accounts for the observed higher-order diffraction in the holography patterns and\nlocal minima near the slit in the self-focusing thermal emission (see the response to Comment 11\nof Review 4 for further details).\nCorresponding Modification:\nWe added the tolerance analysis in Supplementary Text S8.\nWe added the description of the tolerance analysis in the Maintext (page 7, lines 280-281): ‘The\ngroove-fabrication tolerances are within the sculpturing precision of state-of-the-art CNC\ntechniques.’ And lines 288-289: ‘A detailed fabrication-tolerance analysis is provided in\nSupplementary Text S8.’\nSome minor suggestions:\na. I suggest in Fig. 2(a) that only a single unit cell of structure is used to demonstrate the mode profile\nfor the two cases of BM and W-RM. Using multiple periods might be misleading.\nb. In line 247, the authors mentioned an order of magnitude longer of the cohernece length over that in\nAlu’s work. I don’t think this is a fair comparison because the loss associated with the spoof SP in this\nwork is much lower.\nResponse:\na. Thanks for your suggestion. In the initial manuscript, the mode profiles in Fig. 2(a) were\nintended to illustrate the similar period for all three modes. The BM mode is supported by\ndeep-subwavelength grooves and behaves as a dark mode below the light cone. The W-RM\n11\nmode shown in the original submission corresponds to the field profile of a single enlarged\ngroove within the long-period BM structure, effectively folding the band into the radiative\nregime. By tuning the geometry of the deep narrow connecting grooves and the large defect\n(radiative) grooves, we can flexibly tailor both the group velocity and the scattering efficiency\nin our design.\nThe S-RM mode represents an extreme case of W-RM, where the connecting-groove depth is\nreduced to zero. We agree that using mode profiles based on the same single-unit-cell\nconfiguration improves clarity in illustrating the modal properties. Accordingly, we have\nupdated the mode-profile figure (Fig. R7) in the revised manuscript.\nFig. R7 Dispersion diagrams of different spoof SPP modes. The periods of the BM, W-RM\nand S-RM are 0.13λ , 0.70λ and 0.87λ respectively. Groove widths and depths are 0.06λ\n0 0, 0, 0\nand 0.13λ for BM, 0.1λ and 0.11λ for S-RM. For W-RM, the width of the center groove is\n0 0 0\nthe same as that in S-RM and other parameters are the same as those in BM. Mode profiles of\na single unit at frequency ω are shown in the right panel.\n0\nb. The spatial coherence enhancement in our manuscript originates from the group velocity\nenhancement, not the Q-factor or Ohmic loss suppression of the material. In Alu’s paper4, the\nQ factor of their device is larger than 100 but only get ~1.6λ spatial coherence. As a\n0\ncomparison, more than 10λ spatial coherence is achieved with only ~30 Q-factor, as\n0\ndiscussed in the response to Comment 3.\nA key structural innovation of our design is the decoupling between the resonant cavity and\nthe emission surface. Consequently, the coherence of the emitted light no longer relies on a\nband-edge standing wave mode. This design freedom allows us to harness non-resonant\ntraveling surface modes on the emission surface to enhance spatial coherence. From this\nperspective, the comparison with Alù’s work is appropriate and helps readers appreciate the\nadvantages of using traveling surface modes over band-edge resonant modes for spatial\ncoherence enhancement.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of\nthe Nature Communications initiative to facilitate training in peer review and to provide appropriate\nrecognition for Early Career Researchers who co-review manuscripts.\nReviewer #3 (Remarks to the Author):\nI have read the present manuscript on thermal emitters with great interest, and commend the authors\nfor their creativity, but have a number of concerns that I will address below. Some are technical, some\nare more conceptual, with which I will start.\nResponse: Thanks for recognizing our work.\n1. My main concern relates to the way this work is presented: supposedly outperforming, in a direct\ncomparison as shown in Fig. 1b, thermal metasurfaces comprising a patterned emitter. However, this\ncomparison to me is not a fair one: this work, and for example reference 26, rely on thermal point\nsources (in this case the funnel), the radiation of which is then manipulated via a metasurface. This is\nconceptually completely different from thermal metasurfaces where the entire patterned surface is\nemitting. In a sense, the collection of thermal emission from the backside, funneling it through a\nsubwavelength hole, and then manipulating it, is still cascading different optical functions (generation,\ninducing spatial coherence, and then shaping the wavefront sequentially). This is completely different\nfrom true thermal metasurfaces, where the emitting electromagnetic mode contains fluctuating sources\nthat populate it, resulting in a single step emission process of designer functionality. I therefore\nstrongly discourage the authors to make such strong claims and by making a direct comparison with\n“standard” thermal metasurfaces, particularly because these results are interesting, but claiming this\nmuch discredits the work.\nResponse: Thank you for your comment. We agree that our meta-emitter differs fundamentally\nfrom widely reported single-layer thermal metasurfaces, in which the patterned emitting surface\nsimultaneously acts as the source of thermal fluctuations and the wavefront-shaping interface.\nWhile such single-layer architectures offer clear advantages in fabrication simplicity and radiative\npower, they inherently face strict limitations when attempting to simultaneously achieve high\ncoherence and complex wavefront shaping due to their limited design degrees of freedom.\nAs summarized in Fig. 1b, existing single-layer thermal metasurfaces have not yet experimentally\ndemonstrated complex thermal wavefront engineering. Decoupling the functions of coherence\nenhancement and wavefront formation—through a separate resonant cavity and a top scattering\nsurface—provides a practical and experimentally validated route toward complex thermal\nemission, as demonstrated in both our work and that of Yu et al.5\nYu’s method, however, encounters experimental challenges, including limited emission power\nfrom a single subwavelength unit and position-dependent absorptive-structure fabrication. In\ncontrast, our design leverages a bottom-surface mode to enhance emission power and utilizes a\nsingle slit to achieve position-controlled excitation of top-surface waves, enabling experimentally\ndetectable and well-defined wavefront shaping.\nAlthough our architecture appears reminiscent of cascading multiple components—analogous to\ntraditional systems that employ spatial, spectral, and polarization filters together with phase plates\n(Fig. R8a)—the underlying physics is distinctly different. In conventional thermal optical systems,\nthese components operate in the far field, where a large portion of energy is inevitably lost at each\nfiltering stage.\nIn our meta-emitter, all processes occur in the near field, allowing the energy to be efficiently\ncoupled between each component. The thickness of our design is about 1λ , which is similar to\n0\nthat of conventional thermal metasurfaces, as shown in Fig. R8b. Furthermore, similar to standard\nthermal metasurfaces, the bottom surface in our design also hosts fluctuating thermal sources,\n13\npopulated through the same physical mechanism. The key distinction is that this energy is guided\nthrough a central resonant tunnel and subsequently scattered into free space, enabling engineered\ncoherence and tailored wavefronts within a unified thermal-emission framework.\nIn addition, the achromatic wavefront-shaping behavior enabled by the out-of-plane wavevector\ndispersion in the resonant tunnel—discussed in Fig. 3c–d of the Main Text and Supplementary\nText S4—is another unique feature of our design. This property allows a broader spectral range to\ncontribute to the desired wavefront, thereby enhancing energy utilization efficiency.\nFinally, we use the term meta-emitter to emphasize the architectural distinction between our\napproach and traditional single-layer thermal metasurfaces. Although our structure introduces\nadditional complexity to enhance emission performance, comparing our results with those of\nconventional metasurfaces provides a meaningful benchmark for readers and highlights the\nbenefits afforded by the added design flexibility.\nFig. R8 The comparison of the traditional optical system and our meta-emitter in thermal\nwavefront shaping. a. Traditional system relies on spatial, spectral, polarization filter and phase\nplate to get wanted thermal wavefront. b. Meta-emitter integrates all these functions together to\nget the desired thermal wavefront.\nCorresponding Modification:\nWe added the above structural and thermal emission mechanism comparison of our design and\nconventional thermal metasurfaces in the Supplementary Text S1.1.\nWe modified the description of our architecture in the Maintext (page 3, lines 115-126): ‘The\nmeta-emitter incorporates near-field synergistically coupled lossy and lossless\nboundaries, which can be aptly conceptualized as a dual-funnel (DF) model, as\nschematically shown in Fig. 1c. Similar to conventional thermal metasurfaces, the bottom\n‘lossy funnel’ engineers the thermally populated electromagnetic mode and transports it\nto the “lossless funnel” through a coupling resonant “tunnel” while extending their lifetime.\nThese long-lifetime photons correlate the scatterers in the “lossless funnel” with designed\nphase delay, thereby efficiently translating the temporal coherence (W) to the spatial\nt\ncoherence (W ) of tailored emission waves. The additional wavefront design flexibility\nx\ncomes from the increased structural complexity. The thickness of the design can be\nsuppressed to 1λ , which is similar to that of conventional thermal metasurfaces.’\n0\n2. Furthermore, very strong claims are made about creating holograms, but in this case the holograms\nappear to simply be stitched-together lenses as in the first example, given that the holograms comprise\nonly straight lines and the sample appears to be a number of combined lenses even in the photo. Again,\n14\nthis is an interesting demonstration, but this is not a hologram in the traditional sense, where a\ncontinuous phase profile yields a structure a certain distance away from the emitter that can comprise\nmore than only straight lines (two cylindrical lenses next to each other also are not typically\nconsidered a hologram).\nResponse: Thanks for your comment. We agree that the presented holography is not a full two-\ndimensional (2D) holographic demonstration, but rather a combination of two one-dimensional\n(1D) holography operations along orthogonal directions. However, we think there’s a little\nmisunderstanding regarding this aspect. Holography techniques are usually utilized as a measure\nto demonstrate the wavefront manipulation capability. It is exhibited in the 2D pattern for 2D\nphase profile and 1D intensity distribution for 1D phase profile.\nFor each 1D holography design in our work, although the analytical phase expression resembles\nthe superposition of two cylindrical lens profiles, the resulting phase along the surface remains a\ncontinuous, well-defined 1D phase distribution, as detailed in Supplementary Text S7 and Figs.\nS10–S14 (with Fig. S12 replotted in Fig. R9). The use of cylindrical-lens-like phase expression is\nsolely for defining the target intensity distribution at the reconstruction plane.\nIn addition to phase modulation, the scattered-light intensity is also modulated in Designs 4–6,\nfurther supporting the validity of our holography implementation. Therefore, the emission phase\nprofile on the top surface follows standard holographic-design principles and constitutes a 1D\nholography and quasi-2D numerical-pattern reconstruction.\nMoreover, we explicitly discuss the design considerations and challenges of realizing a true 2D\nhologram in the second paragraph of the Discussion section and in Supplementary Text S11.\nFig. R9. Simulated and measured intensity profiles for designs 1-6.\nCorresponding Modification:\nWe clarified the description of our holography demonstrations by explicitly introducing the term\n“quasi-2D holography” to distinguish it from true 2D holographic implementations.\nSome technical concerns:\n3. - One big practical difference is that thermal metasurfaces can reach unity emissivity for specific\npolarization states, but in this case the subwavelength funnel fundamentally limits the amount of power\nthat can be emitted (a single mode waveguide is also bound by a blackbody limit, it cannot transmit\nmore power than that). Thermal metasurfaces are already relatively weak emitters even at unity\n15\nemissivity (compared to coherent sources), so I think this is important to address, especially if you have\nonly one funnel per many wavelengths: how does this affect the emitted power?\nResponse: This is an excellent and critical question for our meta-emitter design, as the emission\npower of thermal metasurfaces is indeed essential for practical applications. Even for structures\nwith unity emissivity, thermal metasurfaces generally remain relatively weak emitters and\ntypically require a large footprint (on the order of 100λ-1000λ) to achieve sufficiently high\nemission power. In our case, the energy efficiency is approximately 10%, as analyzed in our\nresponse to Comment 3 of Reviewer #1. To further enhance the emission power, we proposed a\nphotonic-crystal defect-cavity scheme based on our current design.\nIncreasing the number of funnels per wavelength is certainly an intuitive approach to boosting the\nemission power. However, the key challenge lies in ensuring coherent resonance among the\nfunnels. If the funnels resonate in an uncorrelated manner, the resulting emission is merely an\nintensity superposition of the individual emission profiles, which increases the total power but\ndoes not yield constructive interference or a more refined emission pattern. Achieving phase\ncorrelation across the funnel array requires an additional bottom surface-mode cavity, which could\nenhance the power while also enlarging the spatial coherence on the top surface. This is a\npromising direction that warrants further exploration in future work.\nCorresponding Modification:\nWe have incorporated the above discussion into the last paragraph of Supplementary Text S3.2\n(“Efficiency optimization strategy,” Supplementary page 18, paragraph 2).\n4. - How is the blackbody modeled and taken into account in determining the efficiency and size of the\ncollection area? As far as I can tell this is not mentioned, which makes me think the black body is not\nincorporated directly – however, if so, that would result in a misleading collection size.\nResponse: Thank you for your question. The blackbody (BB) coating acts as a weakly absorptive\nlayer that provides thermal photons in the experiment. It is not explicitly included in the initial\nsimulations; instead, we directly use random dipoles to model thermal excitation. The effective\ncollection size is determined by the decay length of the surface mode, which depends on the loss\nproperties at the bottom surface. Below, we provide additional analysis.\nThe blackbody paint used in our experiment (Musou Black 2.0) is highly absorptive in the visible\nand infrared bands but only weakly absorbing in the millimeter-wave regime. We measured the\nreflectivity of a ~500 μm BB layer coated on a copper plate using a vector network analyzer\n(VNA). The layer thickness was determined by measuring the height difference between the bare\nand the coated metal surfaces with a micrometer. As shown in Fig. R10a (measurement setup in\nthe upper inset), the reflectivity decreases only slightly to ~96% compared with the bare copper\nsurface. To characterize the absorption of this material, we simulated a 500 μm BB layer on metal.\nBy setting the extinction coefficient (imaginary part of refractive index) to 0.02, the simulated\nreflectivity also reduces to ~96%. The corresponding loss density distribution is shown in the inset\nof Fig. R10a. These results confirm that the BB layer behaves as a weak absorber at millimeter-\nwave frequencies.\nWe further repeated the lifetime and decay-length simulations for the three surface modes\ndiscussed in the main text, now including a ~60 μm BB layer (the same thickness used in the\nexperiment). The results are presented in Fig. R10b, with the corresponding loss distributions in\nand W-RM modes by approximately two and one orders of magnitude, respectively, compared\nwith those shown in Fig. 2g. Because the S-RM mode features weak localization, this additional\nloss has only a minor effect—its lifetime and decay length decrease by only ~10% relative to the\nlossless case. Considering the size of our device, only the S-RM mode is used in the experiment,\nand the BB layer therefore has negligible impact. For larger devices, however, such absorption\nwould indeed reduce the effective collection area.\nFig. R10 Degradation of surface-mode lifetime and decay length due to the BB layer. a.\nReflectivity of BB-coated metal surface. The upper inset shows the measurement setup. Left\nbottom panel presents the loss density distribution. The extinction coefficient of BB is 0.02 in the\nsimulation. b. Lifetime and decay length of the three discussed modes with ~60 μm BB coated. c.\nLoss distribution for the three discussed modes.\nCorresponding Modification:\nWe added the description of this BB layer in the page 7, lines 282-288 : ‘A weak absorptive\nlayer, ~60μm blackbody paint (BB paint, n ≈ 1+0.02i @100GHz), was applied to the\nbottom surface for thermal photon generation, as shown in the inset of Fig. 3h. This BB\nlayer has only a minor influence on the lifetime and decay length of the S-RM mode due\nto its weak localization. A detailed analysis of the lifetime and decay length degradation\nfor surface modes caused by the BB layer is presented in the Supplementary Text S7.’\nWe added the detailed analysis shown above in the Supplementary Text S7.\n5.- How does this extend to the near- or mid-infrared? Thermal radiation at 100 GHz is weak, and\nthermal sources are sought after in the mid-IR. This is an interesting demonstration but has little\npractical implication: at high frequencies, one would require a very lossless material combined with a\nsingle point thermal emitter, similar to the aforementioned reference.\nResponse: This is a practical consideration for utilizing thermal radiation in applications. The peak\nradiation for room-temperature objects lies in the mid-infrared (mid-IR) band. Although the\nconductivity of metals is lower in this band compared with the millimeter-wave regime, it is still\nsufficient for the proposed design. We performed a scaled thermal-focusing simulation, presented\nin Supplementary Text S16 and Fig. S35 (reproduced as Fig. R11 here). The simulation uses the\noptical constants of gold in the mid-IR (Fig. R11a).\nThe scaled device preserves single-mode properties across a broad spectral range (2–10 μm), as\nshown in Fig. S35c. The cross-sectional intensity profile at the focal plane is illustrated in Fig.\nS35d, with a full width at half maximum (FWHM) of ~3.6 μm, close to the diffraction limit (~3.1\n17\nμm), indicating strong spatial coherence along the surface. This design principle can also be\nextended to two-dimensional configurations, such as a radially sculpted bull’s-eye structure; an\nexample is provided in Supplementary Text S15.\nFig. R11. Self-focusing thermal emission design at the mid-infrared band. a. Real and\nimaginary parts of the refractive index of gold used in the simulation. b. The eigenmode profile of\nthe designed structure. Scalar bar: 5 μm. c. Normalized transmission of the device with plane\nwave excitation. d. The cross-section of the focal spot. The FWHM is marked in the plot.\n6- What is the spot size for the NA of the lens? It apparently is near-diffraction limited, but the value\nthat one would expect in the ideal case is not mentioned as far as I can tell.\nResponse: In our study, the designed numerical aperture (NA) of the self-focusing device is\napproximately 0.7, corresponding to an ideal diffraction-limited focal spot of ~0.51λ/N.A. ≈ 0.72\nλ, which assumes a uniform amplitude across the aperture. In our design, we optimized the\nfocusing efficiency at the focal spot, which requires scattering all the surface wave energy at the\nedges of the device. This leads to a nonuniform amplitude distribution across the aperture, slightly\nbroadening the focal line in eigenmode simulations to ~0.97λ, representing the optimal trade-off\nfor focusing efficiency. Consequently, the effective NA is ~0.53 for this FWHM. If one prioritizes\nuniform amplitude rather than focusing efficiency, the scattered wavefront can be engineered to be\nnear-uniform, at the expense of discarding surface wave energy at the device edges, approaching\nthe ideal diffraction-limited focal spot of 0.72λ.\nCorresponding Modification:\nIn page 6, paragraph 2, lines 243-250. we modified the description of N.A. and spatial coherence\nof self-focusing device as follows: ‘It achieves a high focusing efficiency of ~95% and a\ntightly confined focal line of 0.93λ (equivalent N.A. = 0.53) realized under the incidence\n0\nof a nonuniform surface wave across the aperture The focal line width can be further\n.\nreduced to the ideal limit of 0.72λ by engineering the radiative loss of the top surface to\n0\nachieve a near-uniform scattered amplitude wavefront, at the expense of dissipating\nsurface wave energy at the edge of the device.’\n7. To conclude, this is an interesting demonstration of an idea very similar to ref. 26, but it is pitched in\na way that is far too strong: it is a compactified set of cascaded elements (thermal source, pinhole,\nwavefront shaper) with very limited demonstrations of wavefront shaping (as far as I can tell only\nfocusing has been shown, including in a demonstration of a “hologram” formed by multiple lenses),\n18\nwhich is not what a thermal metasurface does, and I therefore don’t think it should be pitched as such\nor compared with them.\nResponse: As we addressed in the responses to Comments 1 and 2, although our structure\nresembles a cascade of source, pinhole, and wavefront-shaping elements, the underlying physics is\nfundamentally different from the traditional case. In conventional thermal optical systems,\nspatial filters (including lenses and pinholes), spectral filters, polarizers, and phase plates operate in the\nfar-field, where most of the energy is lost at each filtering stage. In contrast, in our meta-emitter,\nall processes occur in the near-field, allowing the energy to be efficiently coupled between\nelements. The thickness of our design is about 1λ , which is similar to the conventional thermal\n0\nmetasurfaces. At the thermal energy collection plane, the spoof surface mode also contains\nfluctuating thermal sources, populated in the same manner as widely reported thermal\nmetasurfaces. The crucial difference is that, instead of emitting directly into free space, our design\ncouples the energy into a slit and emits it from the opposite side. The fundamental physics remains\nconsistent with thermal metasurfaces.\nAdditionally, the achromatic wavefront-shaping behavior due to out-of-plane wavevector\ndispersion in the resonant tunnel is a novel feature of our design, which significantly enhances the\ntransmission energy across a broad spectral band. As discussed in Fig. 3c–d of the Main Text and\nSupplementary Text S4, this allows more frequency components to contribute to the desired\nwavefront, thereby improving energy utilization efficiency.\nHolography serves as a tool to evaluate wavefront manipulation capability. It manifests as a 2D\npattern for a 2D phase profile or a 1D intensity distribution for a 1D phase profile. In our design,\neach 1D holography case, although the phase profile may resemble the combination of two\ncylindrical lenses, is in fact a continuous phase along the surface, as shown in Supplementary Text\nS7 and Fig. S10–S14 (Fig. S12 is reproduced in Fig. R8). Beyond phase modulation, the intensity\nof scattered light is also carefully designed in Designs 4–6. It is not simply a superposition of two\ncylindrical lens patterns but involves further optimization of the scattered intensity and phase at\neach pixel. Specifically, additional bias phases φ and ψ are applied in the interfered surface wave\nand triple-focusing phase expression, respectively, to produce a smoother combined phase profile\nand achieve the desired energy distribution. Therefore, the emission phase profile on the top\nsurface follows the holography design strategy and can be regarded as 1D holography or quasi-2D\nholography in terms of numerical demonstration.\nThe corresponding modifications have been incorporated in the responses to Comments 1 and 2.\nReviewer #4 (Remarks to the Author):\nThe authors discuss and experimentally demonstrate a bilayer device which performs sophisticated\ntailoring of the wavefront of thermal emission in the millimeter-wave range. The working principle\nrelies on separately engineering two cascaded layers, coupled via a narrow channel, to achieve more\nflexibility in simultaneously collecting thermal emission efficiently and radiate it with a desired\nwavefront.\nThe approach proposed by the authors is interesting, and, to the best of my knowledge, is conceptually\nnovel. I can recommend publication after the authors have addressed the following comments:\nResponse: Thanks for your positive comment.\n19\n1. As a general comment, there are several information, concepts and definitions which are relegated\nto the supplementary material, but which I think should be mentioned in the main text in order to\nprovide a clearer description (note that the SM is 37 pages long). I provide some example in the next\ncomments.\nResponse: Thank you for your constructive comment. In the initial manuscript, many design and\ndefinition details were included in the Supplementary Materials and Methods section to maintain\nconciseness in the main text. The information mentioned below has now been added to the main\nmanuscript.\n2. Can the authors provide a clear definition for the parameters tau_surf and L_D introduced at lines\n141-142?\nResponse: τ is the surface mode lifetime, which is defined by Q/ω, where Q-factor is obtained\nsurf\nin the eigenmode simulation in COMSOL at the designed frequency. It reflects how long the\nenergy of the mode decays to its 1/e. L is the decay length (propagation length) of the surface\nD\n(cid:2914)(cid:2985)\nmode, which is calculated by 𝜏 ∙𝑣 , where 𝑣 (cid:3404) @ω is the group velocity at the working\n(cid:2929)(cid:2931)(cid:2928)(cid:2916) (cid:2917) (cid:2917) 0\n(cid:3031)(cid:3038)\nwavelength, which reflects how long the distance of the energy of the mode decays to its 1/e.\nCorresponding Modification:\nWe put this definition in the caption of Fig. 2b in the original manuscript. We added this part of\ndescriptions into the manuscript (page 5, lines 176-180) as: ‘As shown in Fig. 2b, with realistic\nmaterial loss, the surface mode lifetime  (defined as  = Q /ω ) and decay length\nsurf surf surf 0\nL (calculated as: 𝐿 (cid:3404)𝜏 ∙𝑣 ) of spoof SPPs can vary in two orders of magnitude,\nD (cid:3005) (cid:2929)(cid:2931)(cid:2928)(cid:2916) (cid:2917)\nproviding the flexibility to accommodate devices in different dimensional scales.’\n3. Can the authors clarify what is the difference between the designs BM and W-RM. Based on the\nsketches in the inset in Fig. 2a, they look the same.\nResponse: Thank you for your comment. In the original manuscript, we presented five periods of\nthe BM mode in Fig. 2a. The W-RM mode can be understood as enlarging the center groove\nwithin these five BM periods, effectively folding the band into the light cone. We apologize for\nany confusion caused by the original figures, which were arranged for conciseness but may have\ninadvertently suggested that the modes were more similar than they actually are, without\nproviding an explicit explanation.\nCorresponding Modification:\nWe changed the mode profile of a single unit for BM in the manuscript and added clarification of\nthe difference between the designed three modes (page 5, lines 166-174): ‘Figure 2a plots the\nband diagrams of three representative spoof SPP modes explored in our design: a\nnonradiative bounded mode (BM) below the light cone and two radiative modes: a weakly\nradiative mode (W-RM) and a strongly radiative mode (S-RM). The W-RM can be viewed\nas a perturbed version of the BM, where enlarging one of the grooves within a long\nperiod folds the band into the light cone. The S-RM represents the extreme case of the\nW-RM, obtained by removing all the interconnection grooves within the period. The band\ndiagrams of three modes with light lines are shown in Supplementary Fig. S5 ‘\n20\n4. Can the authors add a light-cone line to Fig. 2a, to better show which modes are dark?\nResponse: Thank you for the comment. Because the periods of the three modes are different, the\nlight-cone lines are different, and it would make the figures complicated if we plotted all these\nlight-cone lines in the figure. We clarified the band position (below or above the light cone) of\neach mode in the manuscript, as shown in the modification in comment 3, and provided the band\ndiagram of each mode with light-cone line in Supplementary Fig. S5 (same as Fig. R12).\nFig. R12 The band diagrams of BM, W-RM and S-RM with light lines. BM is the dark mode,\nwhich is below the light line, and the two radiation modes (W-RM and S-RM) are above the light\nline.\n5. Line 152, “The energy transfer efficiency from source to radiation interface achieves resonant\namplification, showing a significant enhancement at matched frequencies compared to off-resonance\nstates (Fig. 2d).”.\nIt is not clear how such “transfer efficiency” is defined, and what is the “source to radiation interface”.\nCan the authors provide a more detailed description of how the plots in Fig.2d are calculated?\nResponse: Sorry that we only provided the detailed transfer efficiency definition and calculation\nin the Extended Method section. This efficiency represents the mode energy transfer efficiency,\ndefined as the ratio of the energy of the top surface mode transmitted through the slit to that of the\nbottom surface mode. It quantifies the efficiency of energy transfer from the bottom surface mode,\nthrough the slit, to the top surface mode, demonstrating that the designed three components\n(surface modes on both sides and the connecting slit) are well coupled. We have elaborately\ndiscussed this efficiency calculation and strategies to enhance it in response to Comment 2 of\nReviewer 1. The description of energy efficiency has also been clarified in the manuscript, as\nindicated in the modification for that comment.\n6. At line 148, there is a typo in the definition of tau_WG=Q*omega_0. I assume the correct version is\nQ/omega0\nResponse: Sorry for the typo made in the manuscript, and thank you for the correction, we have\ncorrected it in the manuscript.\n7. I am not sure I fully understand the calculations in Figs. 2f-g. Are these calculations done for a\nbilayer device, or for a single layer? If my understanding is correct, here the authors calculate the\neigenmode of the structure (either single layer or bilayer, and for a fixed design type), and then\nevaluate the function W_zz for a fixed value of z=\\lambda0/10. Can the authors confirm if it is a\nbilayer or single-layer simulation, and if it’s an eigenmode calculation?\nIntuitively, I am not sure I understand why the spatial coherence should show such a strong\ndependence on the parameter W1, since W1 affects only the behavior near the center. Can the authors\nelaborate more on this?\ndesign. The bottom layer (funnel) collects thermal photons, where its spatial coherence primarily\ndetermines the area (or length, in the 1D case) over which thermal photons can be harvested. The\ncentral tunnel extends the photon lifetime via Fabry–Pérot resonance. The top surface then\nconverts this temporal coherence into spatial coherence and scatters the photons with the designed\nphase. Therefore, the spatial coherence of the scattered light mainly depends on the slit and the top\nsurface mode.\nAlthough Fig. 2f shows simulations of the full bilayer structure, the cross-spectral density\ncomponent W appears to be influenced primarily by the top layer. This is because eigenmode\nzz\nsimulations capture the field distribution at a single eigenfrequency. To obtain Fig. 2f, we placed a\ndipole on the bottom surface and simulated the field distribution across the detection band (75–\nBecause the spatial coherence (𝐿 (cid:3406)𝜏∙𝑣 ) is dependent on the total lifetime of slit mode and top\n(cid:2917)\nsurface mode (1/τ = 1/τ + 1/τ ) and the dispersion property (group velocity 𝑣 at the\nWG surf (cid:2917)\noperational frequency) of top surface mode. Slit width w directly affects radiation loss of the slit\n1\nmode and controls the lifetime of slit resonance, as we discussed in the Supplementary Text S14,\nwhich not only affect the field distribution around the slit, but also the correlation behavior of the\nscattered light away from the slit.\nFig. 2g in the original manuscript is calculated by substituting the group velocity and lifetime of 3\nsurface modes at the operational frequency, together with the w dependent slit mode lifetime\n1\nfrom Fig. 2c into the above equations.\nCorresponding Modification:\nWe elaborated temporal to spatial coherence mechanism in the manuscript as the modification\nmade in Comment 3 of Reviewer 1.\nWe added more simulation details in the Method Section in the Main Text.\n8. For the calculations shown in Fig. 3b, the authors need to add more detail in the main text, such as:\nmaterials, material parameters, geometrical parameters, and value of lambda0\nResponse: The material in the simulation is set as copper with conductivity of σ = 5.8×107 S/m.\nOperation wavelength is λ ≈ 3 mm (frequency is 100 GHz). The geometrical parameters of self-\n0\nfocusing design are schematically shown in Fig. R13. The specific parameters are listed in the\nTable R1.\nFig. R13 Structure and parameters in 1D self-focusing meta-emitter. The variables t, h and p\ni i, i\ndenote the width, depth, and period of the groove structures, with i = 1, 2, 3 and 4 representing the\nregion of the bottom surface (blue zone), the coupling region (S , orange region), the radiative\nc\nsurface (S, green region), and between the local scatters (S, purple region), respectively. The\nr r\nvariables l and w denote the length and width of the WG, respectively.\nparameter Value (mm) parameter Value (mm)\nt 0.35 p 0.45\n1 2\nh 0.4 t 0.376\n1 3\np 2.88 h 0.483\n1 3\n1.15\nk / k t 0.22\nsspp 0 4\n(dimensionless)\nw 0.1 h 0.308\n4\nl 2.9 dh 0.046\n4\nt 0.2 p 0.52\n2 4\nh 0.4\nWe added the material and material parameters in the caption of Fig. 3b, and added the operational\nwavelength at page 6, lines 241-243: ‘For the optimized structure, the eigenmode field\nanalysis (Fig. 3b) reveals a well-defined focusing pattern with a focal length f = 10λ (λ ≈\n0 0\n3 mm) and numerical aperture N.A. = 0.7.’\nBecause of the word limit on the main text and figure captions, we didn’t list the specific\ngeometric parameters there. These details are provided in the Supplementary Text S4 and Table\nS2.\nWe modified the description of Supplementary Text S4 in the Main Text (page 6, lines 239-241)\nas: ‘Details on parametric optimizations, including k , groove dimensions, and S\nsspp c\nfunctionality, are provided in Supplementary Text S4 and Table S2.’\n9. Line 207, “It achieves a high energy focusing efficiency of ~95%”. Please clarify how focusing\nefficiency is defined here\nResponse: Sorry that this efficiency definition and optimization were presented in the\nSupplementary Text S2 before. This focusing efficiency is defined as the ratio of the integrated\nenergy at the focal line (indicated by the short black line in Fig. S3b) to the energy at the WG exit\n(depicted by the white solid line in the inset of Fig. S3b). High focusing efficiency in our design\nindicates that almost all of the emitted energy contributes to the focal lines. Low background noise\nand sharp focusing properties can be obtained.\nCorresponding Modification:\nWe added the focusing efficiency definition as (page 7, lines 243-247): ‘It achieves a high\nfocusing efficiency of ~95% (defined as the the ratio of the integrated energy at the focal\nline to the energy at the WG exit, details can be found in the Supplementary Text S4 and\nFig. S6) and a tightly confined focal line of 0.93λ\n0\n(equivalent N.A. = 0.53) realized under\nthe incidence of nonuniform surface wave across the aperture\n.’\n10. For the simulations in Fig. 3e, the authors say that “we performed 100 simulations with randomly\ndistributed dipole excitation to obtain statistical results”. Can the authors clarify in which region of\nspace the dipoles are placed? Are the dipoles only placed in the region below the bull’s eye structure?\nIf yes, is this a realistic assumption? In a realistic scenario, any part of the device will emit thermal\nradiation.\nResponse: In the simulation, an electric dipole with a random phase, orientation angle, and x, y\ncoordinates is placed below the bottom surface at a distance of 0.3λ (0.9 mm) as the source to\ncapture the intensity profiles at different frequencies for each plane. This simulation is repeated\n100 times, each with a different random dipole source. The average of all simulated intensity\nprofiles is then calculated to obtain the statistical emission profiles for each plane, as shown in Fig.\n3e-g. Our sample is made by copper with high conductivity in the microwave band, whose thermal\nradiation can be ignored. The thermal radiation is mainly caused by the blackbody paint under the\nbottom surface, which can be simulated by placing dipoles below the bottom surface in the\nsimulation. The parameters of the black body and its absorptivity are analyzed in the response to\nComment 4 of Reviewer 3.\nCorresponding Modification:\nWe added the dipole position information in the manuscript (page 7, lines 263-266): ‘To\nnumerically characterize the stochastic thermal emission behavior of our design, we\nperformed 100 simulations with randomly distributed dipole excitation below the bottom\nsurface at a distance of 0.3λ to obtain statistical results (detailed simulation methods\n0\ngiven in Method).’\n11. I notice some difference between experiments (Fig. 3i) and simulations (Fig. 3e). In particular, in\nthe experimental data there is an intensity local minima near the center of the structure (close to the\nWG element), while the simulations predict a high intensity. Can the authors comment on this?\nResponse: We found that the intensity of local minima near the center slit originates from the\nnon-ideal interference of waves scattered from the left and right sides of the lens grooves, caused\nby depth fabrication errors in the coupling region around the slit exit. We conducted a simulation\nwith a 10% larger groove depth; in this case, the scattered waves from the first grooves on both\nsides of the device can form destructive interference in the near field, as shown in Fig. R14b. Fig.\nR14a (same as Fig. 3e) represents the ideal scattering case, where the two grooves constructively\ninterfere at the center. The measured intensity profile in Fig. R14c resembles a combination of\nboth cases, with the interference minima more pronounced in the near field. This fabrication error\nmainly affects the near-field intensity distribution and has a minor impact on the far-field pattern,\nconsistent with our experimental observations.\nFig. R14 Constructive and destructive interference cases around the slit exit. a-b. Simulation\nresults for the ideal structural parameters and non-ideal structural parameters with 10% deeper\ngroove depth of the coupling region (S in Fig. 3a in the Main Text). c. Measured intensity profile.\nc\nCorresponding Modification:\nWe added the above discussion in the Supplementary Text S8.4.\nWe added the description of this difference on page 8, lines 301-305: ‘The intensity local\nminima near the center WG originates from the non-ideal interference of the waves\n24\nscattered from the left and right lens grooves, due to the depth fabrication error in the\ncoupling region around the slit exit. The detailed analysis can be found in Supplementary\nText S8.4.’\n12. When describing the fabricated structure in the caption of Fig. 3h, the authors mention that the\nbottom side is “coated with a blackbody layer”. The need for this layer is not mentioned elsewhere in\nthe paper, and it’s not clear which role it plays. Is the presence of this layer captured in the simulations?\nPlease include a more detailed description of it\nResponse: Sorry that we didn’t explain the function of this layer directly in the original\nmanuscript. The blackbody coating serves as a weak absorptive layer that provides thermal\nphotons in the experiment, which corresponds to the bottom structures in the model in Fig. 1a-b.\nWe made a detailed analysis of this layer, and the corresponding modifications are provided in the\nresponse to Comment 4 of Reviewer 3.\nSincerely Yours,\nOn behalf of all authors,\nYungui Ma, Professor, College of Optical Science and Engineering, Zhejiang University,\nChina\nCheng-Wei Qiu, Professor, Department of Electrical and Computer Engineering,\nNational University of Singapore, Singapore\nReference:\n1. Janssen, O. T. A., Urbach, H. P. & ’T Hooft, G. W. Giant Optical Transmission of a\nSubwavelength Slit Optimized Using the Magnetic Field Phase. Phys. Rev. Lett. 99, 043902\n(2007).\n2. Fan, S., Suh, W. & Joannopoulos, J. D. Temporal coupled-mode theory for the Fano\nresonance in optical resonators. J. Opt. Soc. Am. A 20, 569 (2003).\n3. Wonjoo Suh, Zheng Wang, & Shanhui Fan. Temporal coupled-mode theory and the\npresence of non-orthogonal modes in lossless multimode cavities. IEEE J. Quantum Electron.\n40, 1511–1518 (2004).\n4. Nolen, J. R., Overvig, A. C., Cotrufo, M. & Alù, A. Local control of polarization and\ngeometric phase in thermal metasurfaces. Nat. Nanotechnol. 19, 1627–1634 (2024).\n5. Zhou, M. et al. Self-Focused Thermal Emission and Holography Realized by\nMesoscopic Thermal Emitters. ACS Photonics 8, 497–504 (2021).\n25"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have reviewed the authors' responses to my previous comments and the corresponding revisions. All concerns have been adequately addressed. I have one optional suggestion for the Discussion section: The strategy for arbitrary control of the thermal emission wavefront fundamentally relies on spoof surface plasmon polariton excitation, which explains the millimeter-wave demonstration. If this design is extended to the mid-infrared regime—where thermal emission peaks according to Planck's law—how would the increased metallic losses affect the proposed capabilities?\nOtherwise, I recommend publication in Nature Communications.", "reviewer_2": "The authors have adequately addressed my concerns; therefore, I recommend accepting this paper.", "reviewer_3": "I appreciate the effort the authors have put into the manuscript in response to my and the other reviewers' comments. However, regarding my main concern, I do not feel this has been addressed entirely. The authors claim that the physics between this platform and a cascaded source/pinhole/metasurface stack is completely different, but I disagree. The role of the tunnel/pinhole is clearly crucial, because it ensures a point source/spatial coherence – and with spatial coherence it is trivial to manipulate wavefronts, which is the backbone of virtually all metasurfaces. As such, in my opinion, this work has not addressed the inherent challenges to thermal metasurfaces, which the authors outline nicely, but instead they have compressed the standard approach (lamp/pinhole/metasurface) into a much smaller footprint.\nDon't get me wrong: I like this work, I think it is impressive and should be published, but it currently is minimizing/dismissing the achievements of the community starting with Greffet who are pursuing control over thermal radiation sculpting by suggesting that this approach is superior. It is not, it is merely a different and smart approach to ensure spatial coherence which enables easier manipulation of wavefronts.\nI would strongly encourage the authors to remove direct comparisons between the two approaches on the basis of e.g. spatial coherence length, since, again, in the present metasurface the coherence is imprinted by the pinhole. Likewise, I would ask the authors to explicitly mention in the last paragraph of the introduction and in the architecture section that the key ingredient in enabling wavefront manipulation is the thermal point source, as in the other reference that proposes thermal holography.\nNote that there are a number of typos in Fig. 1 and the caption.", "reviewer_4": "The authors have addressed all my comments satisfactorily. I can now recommend this manuscript for publication in Nature Communications."}, "author_response": "To Reviewers:\nReviewer #1:\nComments:\nI have reviewed the authors' responses to my previous comments and the corresponding revisions. All\nconcerns have been adequately addressed. I have one optional suggestion for the Discussion section:\nThe strategy for arbitrary control of the thermal emission wavefront fundamentally relies on spoof\nsurface plasmon polariton excitation, which explains the millimeter-wave demonstration. If this design\nis extended to the mid-infrared regime—where thermal emission peaks according to Planck's law—\nhow would the increased metallic losses affect the proposed capabilities?\nOtherwise, I recommend publication in Nature Communications.\nResponse: Thank you for your acknowledgment and recommendation. This is a good question\nregarding the extension of our design to the widely used mid-infrared (mid-IR) regime, where\nthermal emission peaks according to Planck’s law.\nIn the mid-IR, increased metallic loss indeed leads to higher absorption for both the spoof surface\nplasmon polariton (SSPP) mode and the connecting slit waveguide. In the demonstrated\nmillimeter-wave regime, the conductivity of copper is taken as σ = 5.8 × 10⁷ S/m at a wavelength\nof 3 mm. By comparison, in the mid-IR, gold exhibits complex refractive index n = n + ik = 2 +\n39i at λ = 5.8 μm, corresponding to an effective conductivity of σ = 2nkωε ≈ 4.6×105 S/m.\nconfinement configuration, such as the S-RM mode discussed in the main text. We simulated a\nstructure scaled by a factor of 600 operating in the mid-IR. In this case, the Q factor decreases\nfrom 98 to 81, corresponding to only a 17% degradation in temporal coherence. The group\nvelocity remains high (≈ 0.93c), comparable to that of the original design, while the spatial\ncoherence degradation is also limited to approximately 17%. Furthermore, the temporal decay of\nthe scaled mode can be compensated by further reducing the size of the radiative grooves\n(perturbations), thereby maintaining a low radiative decay rate for the surface mode.\nFor the slit waveguide, the Q factor decreases from 31 to 21, which primarily affects the spatial\ncoherence of emission from the top surface. This radiative loss can be substantially suppressed\nthrough geometric optimization of the slit, for example by enlarging the chamber region to reduce\nelectromagnetic energy density near the metal walls. This mechanism is discussed in detail in\nSupplementary Text S14 and Fig. S28.\nFinally, a fully scaled model operating at 5.7 μm is simulated in Supplementary Text S16. The\nresults demonstrate single-mode operation over a broad wavelength range (2–10 μm) and near–\ndiffraction-limited focusing performance, with a full width at half maximum (FWHM) of\napproximately 3.6 μm.\nCorresponding Modification:\n1. In the Discussion part, lines 392-394, we added a brief discussion of the effect of increased\nmetallic absorption loss on the thermal meta-emitter design as: ‘The increased absorption\nloss can be mitigated by employing low-confinement surface mode, such as the\nS-RM mode in Fig. 2a, and compensated by further reducing the size of the\nradiative perturbations.’\n1\n2. We have added the above discussion into the Supplementary Text S16 (Supplementary pages\n48-49, paragraphs 1-2).\nReviewer #3:\nComments:\nI appreciate the effort the authors have put into the manuscript in response to my and the other\nreviewers' comments. However, regarding my main concern, I do not feel this has been addressed\nentirely. The authors claim that the physics between this platform and a cascaded\nsource/pinhole/metasurface stack is completely different, but I disagree. The role of the tunnel/pinhole\nis clearly crucial, because it ensures a point source/spatial coherence – and with spatial coherence it is\ntrivial to manipulate wavefronts, which is the backbone of virtually all metasurfaces. As such, in my\nopinion, this work has not addressed the inherent challenges to thermal metasurfaces, which the\nauthors outline nicely, but instead they have compressed the standard approach\n(lamp/pinhole/metasurface) into a much smaller footprint.\nDon't get me wrong: I like this work, I think it is impressive and should be published, but it currently is\nminimizing/dismissing the achievements of the community starting with Greffet who are pursuing\ncontrol over thermal radiation sculpting by suggesting that this approach is superior. It is not, it is\nmerely a different and smart approach to ensure spatial coherence which enables easier manipulation\nof wavefronts.\nI would strongly encourage the authors to remove direct comparisons between the two approaches on\nthe basis of e.g. spatial coherence length, since, again, in the present metasurface the coherence is\nimprinted by the pinhole. Likewise, I would ask the authors to explicitly mention in the last paragraph\nof the introduction and in the architecture section that the key ingredient in enabling wavefront\nmanipulation is the thermal point source, as in the other reference that proposes thermal holography.\nNote that there are a number of typos in Fig. 1 and the caption.\nResponse: Thanks for your acknowledgment and suggestions. We didn’t intend to diminish the\nsignificant achievements of traditional metasurfaces that directly sculpt thermal radiation with\nsingle-layer structures. To avoid any potential misunderstanding, we agree to remove this part of\ncomparison and to elaborate more explicitly on the role of the central slit.\nCorresponding Modification:\n1. We have removed the comparison table (Fig. 1b in the original manuscript) and the\ncorresponding comparison descriptions.\n2. In last paragraph of Introduction, page 3, paragraph 2, lines 99-101, we added the explicitly\ndescription of connecting waveguide as:\n“The single connecting waveguide ensures a fixed scattering-phase relationship on the\noutput surface, thereby endowing meta-emitters with arbitrary phase modulation\ncapability.”\nBesides, in the first paragraph of ‘Architecture of Meta-Emitter’ section (page 4, paragraph 1),\nlines 125-127, we added the following clarification:\n“The connecting tunnel serves as an effective coherent source, fixing the scattering\nphase on the lossless funnel, similar to the point source in Ref. 26.”\n2\n3. We checked and corrected the typos in Fig. 1 and its caption.\nSincerely Yours,\nOn behalf of all authors,\nYungui Ma, Professor, College of Optical Science and Engineering, Zhejiang University,\nChina\nCheng-Wei Qiu, Professor, Department of Electrical and Computer Engineering,\nNational University of Singapore, Singapore\n3"}]

Decision letter