All-optical hybrid metasurfaces for ultrafast computational spectrometer and single-pixel imaging

Domain: physics
Published: 2025-12-14
Editorial outcome: minor_revision

Abstract
Terahertz(THz)spectrumholdsimmensepotentialforapplicationsinchemicalsensing,biomedicaldiagnostics,andsecurityscreening.However,conventionalTHzspectrometersrelyonmechanicaldelaylinesandcoherentdetection,posingmajorchallengesforminiaturization,speed,andcost.Here,wereportanactivehybridmetasurfacethatfunctionsasalow-correlationandultrafastspectral/spatiallightmodulatorforTHzcomputationalspectrometerandsingle-pixelimaging.Dispersion-engineeredmetasurfacewithmultiplequasi-boundstatesinthecontinuumachievesdensehighqualityfactorresonancesacrossaspectralrange(0.30-0.55THz),enablingdynamicmodulationthroughall-opticalexcitationofintegratedsiliconpatches.Temporallymodulatedmetasurfacesformameasurementmatrixforincoherentsingle-pixeldetection,allowingnanosecond-scalespectralreconstructionwitharesolutionof0.03THzcoveringa0.25THzbandwidth.Extendingthisplatformtoapixelatedmetasurfacearray,wedemonstrateultrafastTHzsingle-pixelimagingwith3×3pixels.ThisCMOS-compatible,circuit-freestrategyenablesultrafast,compact,andhigh-resolutionTHzspectroscopyandimaging,andofferascalablepathwaytowardnext-generationTHzphotonicdevices.

Introduction
Terahertz (THz) waves hold significant promise for applications in chemical analysis, biological sensing, medical diagnostics, and security screening due to their unique properties including strong penetrability, nonionizing photon energy, and rich spectral fingerprints spanningawiderangeofmolecularvibrationsandresonances1–4.THzspectrometersandimagingsystemsarethusindispensabletoolsforspectralcharacterizationandnondestructiveinspection.WhileconventionalTHztime-domainspectrometers(THz-TDS)enablethedirectacquisitionofbroadbandspectrawithphaseinformation,thesesystemstypicallyrelyoncoherentdetectionschemes(Fig.1a),whichnecessitatebulkyopticalsetupsincorporatingmechanicaldelaylines5–8.Moreover,spatialinformationmustbeacquiredthroughrasterscanning,especiallywhensingle-pixeldetectors(SPDs)areusedforimaging.Asaresult,therealizationofcompact,low-costandportableTHzspectrometersremainsalongstandingchallenge.

ApromisingroutetominiaturizationinvolvesreplacingcoherentdetectionwithanincoherentSPD(Fig.1b),sacrificingdirectaccesstobroadbandphase-resolvedspectraldata.Tocompensate,computationaltechniquessuchascompressedsensing,regularization,anddeeplearningalgorithmsareemployedtoreconstructspectralandspatialinformationfromintensity-onlymeasurements9–11.ThesemethodsrelyonsolvingtheinverseproblemI=ΦX,wherеФdenotesthemeasurementmatrix,Xrepresentstheobject(spectralorspatialdistribution),andIistherecordedintensityvector.Crucially,successfulreconstructiondependsonthemeasurementmatrixcomprisingsufficientlyuncorrelatedspectralorspatialencodingpatterns.

Inthiscontext,spatialorspectrallightmodulators(SLMs)playacentralrolebygeneratingdecorrelatedmodulationframesthatdefinethemeasurementmatrix12–14.Thedegreeofframe-to-framedecorrelationdirectlyaffectsreconstructionfidelity,alongsidefactorssuchasmodulationspeedandspectralbandwidththatdetermineoverallperformance.Bothnaturalmaterials(e.g.,quantumdot,liquidcrystals,andperovskites)andartificialphotonicmaterials(e.g.,metasurfacesandphotoniccrystals)havebeenexploredascandidatesforencodingpatterns,particularlyinthevisibleandinfraredspectralregimes15–19.However,theinherentlyweakinteractionbetweenTHzradiationandmostnaturalmaterialsposesasubstantialbarriertodevelopingTHz-compatibleSLMsthatcombinelowcorrelationcoefficients,broadbandwidth,andhigh-speedmodulation.

Metasurfaceshaveemergedasapowerfulplatformforenhancinglight-matterinteractionsinTHzSLMs20–22.Asartificialmaterials,metasurfacesofferunparalleledflexibilityintailoringthespatialdistributionsandshapingspectralresponsesofelectromagneticfields,enablinglocalfieldconfinementinthespatialdomainandenhancedinteractioninthetemporaldomain23–25.Whenintegratedwithfunctionalnaturalmaterials,suchassemi-conductorsorphase-changemedia,hybridmetasurfacesexhibitdynamictunabilityandserveasversatileplatformsforreconfigurablemetadevicesinmoldingwavefrontandexploringphysicalphenomena26,27.

Inconventionalactivemetasurfaces,enhancedlight-matterinteractionsarenaturallyassociatedwithhighqualityfactorresonances(Q),whereQ≈f/Δf,withfdenotingthecentralfrequencyandΔfthefullwidthathalfmaximum(FWHM)linewidth,andthuspossessaverynarrowbandwidth.High-Qresonancesareoftenleveragedtoreducemodulationthresholdsandamplifytunability.However,theirnarrowspectralwidthinherentlylimitsthebroadbandspectralencodingneededforcomputationalspectrometry.Furthermore,themodulationspeedoftraditionalmetasurface-basedSLMsisconstrainedtothemillisecondrangebytheintrinsicrelaxationtimesofconstituentmaterials,suchasliquidcrystals,phase-changematerials,andmicro-electromechanicalsystems(MEMS),aswellastheRCtimeconstantsofdrivecircuitry28–31.All-opticalmodulationschemesprovideanalternativeroutebyeliminatingcomplexelectronics,therebyenablingultrafastmodulationonthepicosecondtimescale32,33.Nevertheless,theinherenttrade-offbetweenlowspectralcorrelation(enabledbyhigh-Qresonances)andbroadbandspectralcoverageremainsafundamentallimitationforbroadbandcomputationalspectroscopy.

Recentadvancesinthephysicsofboundstatesinthecontinuum(BICs)offerasystematicframeworkforovercomingthistrade-off34–36.BICssupporttheoreticallyinfiniteQfactorsandallowfortheengineeringofradiationchannels,dispersionrelations,andmodalprofiles37–39.PixelatedBICsmetasurfaceshavebeenshowntomaintainhigh-Qresonanceswhileexpandingthespectralcoverage,provingbeneficialforapplicationssuchasultrasensitivehyperspectralimagingandmolecularbarcoding40,41.However,insuchpixelatedschemeswithNspatialelements,eachmetasurfacepixelinteractswithonly1/Noftheincidentpoweratitsresonantfrequency,withtheremaining(N−1)/Nconstitutingbackgroundenergy.Thisspatialmultiplexingapproachinevitablyreducesthesignal-to-noiseratio(SNR).

Toaddressthislimitation,weintroduceastrategythatfoldsmultipleBICsintoabroadbandspectralregion,therebypreservinghighQ-factorswhileenhancinglightutilizationwithoutreducingSNR42.BybreakingsymmetrytoconvertBICsintoquasi-BICs(q-BICs),adensesetofhigh-Qresonancescanbeexcitedacrossawidefrequencyrange.Theseresonancesexhibitdiversemodeprofilesandenhancedsensitivity,enablingbroadbandspectralmodulationwithintrinsicallylowcorrelationcoefficients.

Here,weexperimentallydemonstrateabroadband,low-correlation,andultrafastSLMbasedonhybridmetasurfacessupportingmultipleq-BICs,andshowtheirapplicationsinTHzcomputationalspectrometer(THz-CS)andTHzsingle-pixelimaging(THz-SPI).Asaproof-of-conceptdemonstration,thebroadbandspectralcorrelationisreducedbytransferringthreeq-BICstotheinterestedspectralrange(0.30-0.55THz),andthemetasurfacesaredesignedtomodulatetransmissionspectrainanall-opticalmannerwithnanosecondscaleresponsetime.IncorporatedintoaTHzsystemwithanincoherentSPDwithoutspectralresolutioncapability,theactivehybridmetasurfacesactingasanSLMenablesultrafastTHz-CSwith0.03THzspectralresolutioncoveringa0.25THzbandwidthand1.9nsreconstructiontime.Theall-opticalhybridmetasurfacearchitectureisreadilyextendableforultrafastsingle-pixelimagingviapixelatedmetasurfacearrays.

Results
Ourdesignemploysmetasurfacescomposedofclassicaldouble-gapsplit-ringresonators(DSRRs),chosenfortheirwell-understoodelectromagneticbehavior.Whenthesymmetryoftheunitcellsisbroken,Fanoresonancesemergeduetointerferencebetweenq-BICsresonancesandthebroadradiativebackground.TheseFanoresonances,knownfortheirsharplineshapesandstrongsensitivity,areparticularlyattractiveforuseinsensingandactivemodulation43,44.Toachievebroadbandoperationwithlowspectralcorrelation,weexcitemultipleFanoresonanceswithdistinctmodesymmetries.ThisisaccomplishedbyengineeringthemodedispersionofDSRRsarrangedinasquarelattice.Whenthelatticeperiodicityisdoubled,twoguidedmodesattheBrillouinzoneedgearefoldedtotheΓpoint(insetinFig.2a).Uponsymmetrybreaking,thesemodesbecomeradiative,producingtwoadditionalq-BICs(labeledq-BICIandq-BICIII)inthevicinityoftheoriginalFanoresonance(q-BICII).ThesenewresonancesretaintheBICcharacter,exhibitinganinversequadraticdependenceofQonthein-plane wavevector45. Further details of the band-folding design and simulations are provided in Supporting Information S3.

We introduce symmetry breaking by laterally displacing the capacitive gap in one of the two DSRRs along the x axis (d=27μm, Fig. 2a inset), which breaks both the C2 symmetry of resonator and lattice symmetry, converting all three BICs into radiative q-BICs. Experimentally, we observe three distinct resonances centered at 0.35THz, 0.44THz, and 0.48THz, with Q-factors of 16.1, 10.5, and 22.5, respectively, in good agreement with simulations. The intrinsic ohmic losses of metallic resonators restrict the increase of Q-factors to a higher value comparable to dielectric metasurfaces. Specifically, the three resonances deliver around 60% relative bandwidth centered at 0.42THz46. Mode profiles, visualized via surface current and local electric field distributions (Fig. 2b), reveal that q-BIC I and II correspond to antiphase and in-phase oscillating current loops in the neighboring two DSRRs, deemed as inductive-capacitive (LC) resonances (anti-LC and parallel-LC) with the localized electric field predominantly concentrated in the capacitive gaps. q-BIC III manifests as a dipole-like mode, with electric fields primarily concentrated at the metallic arm edges.

These distinct local field distributions underpin our design of hybrid metasurfaces for broadband low-correlation modulation. To achieve active modulation, we embed silicon patches at strategic locations within the resonators. Optical pumping of the silicon alters its conductivity with excess photocarrier accumulation, thereby modifying the equivalent inductance and capacitance on nanosecond timescales47,48, which significantly impacts the resonance conditions due to strong local field enhancement. Detailed influence on each q-BICs resonances due to perturbation of silicon conductivity is summarized in Supporting Information S4, and perturbation in capacitive gaps tunes q-BIC I and q-BIC II, while perturbation in resonator arms mainly contributes to the modulation of q-BIC III. Hybrid metallic split-ring resonators (SRRs) integrated with semiconductors have been studied for various applications in the THz regime49. Leveraging their well-established physics, we use band folding to broaden the modal spectrum and apply it to THz-CS/SPI. More importantly, our theoretical framework based on q-BICs and dispersion engineering is general and applies to a wider range of metallic or dielectric resonators, rather than being limited to DSRRs.

We investigate two hybrid metasurface types to assess their modulation correlation: frequency-shift type, where silicon patches are embedded in the capacitive gaps without electrically shorting them, thus preserving mode profiles while shifting resonance frequencies; and mode-change type, where gaps are fully filled in one resonator, leading to quenching of LC resonances upon conductivity increase. Spectral modulation characteristics are presented in Fig. 2c. In the frequency-shift design, all three q-BICs exhibit redshifts and reduced depths with increasing silicon conductivity. In contrast, the mode-change design suppresses the original resonances and induces new dipolar responses. Correlation coefficients between modulated and unmodulated spectra reveal that mode-change designs consistently produce lower values, indicative of greater spectral diversity and lower inter-frame redundancy (Fig. 2d).

Having addressed the trade-off between broadband spectral coverage and high-Q resonances to reduce modulation correlation coefficients, we next turn to the implementation of an ultrafast measurement matrix essential for high-speed computational reconstruction. To this end, we fabricated the metallic DSRRs and integrated them with a patterned silicon epilayer on a sapphire substrate. Optical micrographs of the fabricated hybrid metasurfaces are shown as insets in Fig. 3a and e (see “Methods” for detailed fabrication procedures). To probe ultrafast broadband spectral modulation, we employed an optical-pump THz-probe (OPTP) system using femtosecond optical pulses with a fluence of 0.52mJ/cm2 offering external pump perturbation. Upon excitation, photocarriers are generated in the silicon patches, modulating the local conductivity and consequently altering the resonant behavior of the metasurface at specific frequencies. We recorded the electric field variation ΔE/E0 at the peak of the time-domain signal pulse, where E0 represents the reference field amplitude without pump measured under identical conditions, and ΔE=Epump - E0. The subsequent excitation and relaxation dynamics of these nonequilibrium carriers were tracked using the delayed THz probe. The photocarrier density reaches its peak ~27ps after the pump pulse (defined as τ = 0ps), and then gradually decays over ~1.9ns via carrier recombination and trapping processes. From the measured differential transmission spectra, we estimate the time-dependent silicon conductivity to vary from ~4800S/m after excitation, decaying to ~300S/m over the 1.9ns relaxation period (see Supporting Information S5).

Although both types of hybrid metasurfaces (frequency-shift and mode-change) exhibit similar temporal carrier dynamics, governed primarily by the intrinsic properties of the silicon epilayer, their spectral evolution diverges significantly due to differences in silicon distribution and its interaction with the local electromagnetic fields. Figure 3b and c illustrate the spectral evolution of the frequency-shift metasurface based on both experimental data and numerical simulations. As photocarriers accumulate, all three q-BIC resonances undergo a redshift and exhibit reduced resonance depths, attributable to the modulation of equivalent inductance, capacitance, and resistance in the resonator circuit. These spectral features recover progressively to their steady-state positions over the 1.9ns timescale. A persistent, unperturbed feature at ~0.41THz, observed experimentally, is attributed to Fabry-Pérot interference within the sapphire substrate, which is relatively insensitive to free-carrier modulation. In contrast, the mode-change metasurface demonstrates a more complex spectral evolution (Fig. 3f and g). At τ = 0ps, a new resonance emerges at ~0.38THz, replacing the original q-BIC I and q-BIC II resonances, suggesting a mode transition induced by gap quenching. The q-BIC III resonance, however, shows a redshift trend similar to that observed in the frequency-shift design. All three resonances gradually return to their unperturbed states over the course of 1.9ns. This ultrafast and dynamic spectral evolution effectively forms a time-varying modulation mask as a measurement matrix enabling object encoding and spectral information reconstruction on a nanosecond timescale.

The presence of three distinct q-BICs resonances induces strong spectral undulations across the broadband range of 0.30 to 0.55THz, providing a robust encoding basis for a THz-CS. Upon optical excitation, the hybrid metasurface generates a sequence of temporally varying spectral frames with high contrast, each exhibiting unique modulation patterns, within a nanosecond timescale. These dynamically evolving frames effectively constitute a measurement matrix for spectral reconstruction. To quantitatively evaluate the encoding quality of the measurement matrix formed by the hybrid metasurfaces, we calculate the spectral auto-correlation function50–52:

C(Δf) = <T(f,τ)T(f+Δf,τ)>f/τ / (<T(f,τ)>f/τ <T(f+Δf,τ)>f/τ)

where T(f,τ) is the transmission amplitude at frequency f and at time τ, and <·>f/τ indicates averaging over frequency f or time τ. The auto-correlation function is normalized such that C(0)=1. The resulting auto-correlation profiles, shown in Fig. 3d and h, characterize the degree of spectral diversity encoded across modulation frames.

The half-width at half-maximum (HWHM) of C(Δf), i.e., the frequency offset where the auto-correlation falls to C(Δf)=0.5, defines the achievable spectral resolution when the active hybrid metasurface is used as an SLM in THz-CS. Based on 10 uniformly spaced temporal frames selected from the 1.9ns carrier relaxation window, we extract spectral resolutions of 0.028THz and 0.024THz for the frequency-shift and mode-change metasurfaces, respectively (see Supporting Information S6). The improved resolution of the mode-change design is attributed to its greater spectral decorrelation across frames, and it is thus selected as the prototype for experimental validation of the THz-CS using an SPD in a simplified system configuration (More modulation strategies are provided in Supporting Information S7).

The proposed simplified THz-CS architecture is illustrated in Fig. 4a. A broadband THz pulse generated from an antenna is transmitted through the sample, forming an unknown spectral response S(f). Unlike conventional systems that rely on mechanical delay lines and coherent detection schemes, the proposed system leverages an active hybrid metasurface as an SLM in conjunction with an SPD, and divides partial femtosecond pulses from transmitter path to drive the SLM, thereby significantly reducing the complexity, size, and cost of the spectrometer.

Spectral modulation frames are generated by photoexciting the hybrid metasurface with a pulsed optical pump, synchronized with the SPD to trigger the series of signal acquisition on a nanosecond timescale. The unknown spectrum is encoded by a measurement matrix T(f), determined by the metasurface’s dynamic response across modulation frames. The SPD records the resulting intensity vector I, from which the original broadband spectrum can be reconstructed computationally. The fidelity of the reconstruction both in spectral resolution and SNR is governed by the spectral resolution and inter-frame correlation of T(f), the number of modulation frames, algorithmic efficiency, and system stability (see Supporting Information S8 and S9). In addition, we analyze the impact of system repeatability and experimental uncertainties on the performance (see Supporting Information S11 and S12).

As a proof of concept, we evaluated the spectral reconstruction performance using a set of THz band-pass filters with different center frequencies and approximately Gaussian passbands. This choice reflects the intended application of our THz-CS, detecting absorption features of biomacromolecules, pharmaceuticals, and gases, whose lineshapes are typically well approximated by Gaussian under practical conditions. The attempts to reconstruct rectangular and triangular waveforms can be found in Supporting Information S13. With only 10 modulation frames, the reconstructed spectra faithfully reproduce the filters’ original lineshapes and central frequencies, achieving a FWHM of ~0.03THz (Fig. 4b), comparable to the resolution obtainable from conventional time-domain spectroscopy with a scan time of ~30ps. The average peak frequency difference (Δf) between reconstructed and original spectra is 3.95×10-3THz, with a minimum of 4×10-5THz. The average frequency resolving power (R = f/Δf) at the specified operating frequency is ~1193, which is comparable with previous reports (Fig. 4c)53,54. Spectral reconstruction with different features was carried out for a single-pass filter with a broader bandwidth and spectrum with two closely spaced Lorentzian peaks (~0.05THz apart), and the reconstructed spectra reveal excellent quality via the proposed THz-CS system (Fig. 4d). We note that the reported 0.25THz operational bandwidth refer to the frequency range over which the resonant low-spectral-correlation modulation is significant and the reconstructed spectra remain accurate. To achieve a broader operational bandwidth, the DSRR metasurface can be geometrically scaled to produce multiple modules, each providing a substantial frequency coverage in distinct spectral regions. By assembling these modules, the system can achieve an ultra-wide operational bandwidth while maintaining low spectral correlation.

The overall quality of the reconstructed spectrum is characterized by the peak signal-to-noise ratio (PSNR, PSNR=10 log10(MAX2/MSE) where MSE is mean squared error, MAX is maximum possible pixel value of the spectrum), and a larger PSNR indicates a higher-quality reconstructed spectrum55. Based on the example of single-peak bandpass filter spectrum centered at 0.39THz, we observe an improved PSNR with a larger number of measurements under both types of hybrid metasurface SLMs (Fig. 4e). The PSNR increases as the number of measurements rises from 5 to 15, but it saturates once the frame count exceeds ~20, indicating that additional frames contribute redundant information due to non-negligible correlation between adjacent modulation states. With 10 measurements, the mode-change SLM attains an excellent reconstruction PSNR of 30.1dB, commonly regarded as the threshold for high-quality reconstruction15,56,57. The comparison between the two types of hybrid metasurface SLMs indicates an overall better performance of mode-change type with lower number of measurements. We also investigate the relationship between the PSNR of the reconstructed spectra and the spectral resolution (see Supporting Information S10).

In addition to spectral reconstruction by using the broadband, low-correlation ultrafast modulation, the hybrid metasurface can be spatially pixelated to encode spatial information, functioning as a high-speed SLM for THz-SPI. In a conventional THz-SPI system, a continuous-wave THz beam transmits through or reflects from the target object X, and is subsequently modulated by predefined spatial masks Φ before being detected by an SPD, forming an intensity vector I (Fig. 5a). However, the overall imaging speed is severely constrained by the slow modulation rate of conventional SLMs, which often rely on electronically driven spatial masks with complex peripheral circuitry. We present a circuit-free, all-optical solution based on our hybrid metasurface, which enables ultrafast spatial modulation by harnessing the rich spectral responses and mode distributions of band-folded q-BICs.

Operating at a central frequency of 0.45THz, we design a 3×3 pixelated mode-change type metasurface array, where each pixel is scaled by a geometrical factor S (Fig. 5b) to induce distinct temporal transmission dynamics in response to a uniform optical pump. The pixels are selected to produce strongly contrasting transient responses during the relaxation process, as detailed in Supporting Information S15. Figure 5c illustrates the time-resolved transmission responses of the nine pixels over a 1.9ns interval, revealing diverse and uncorrelated temporal modulation patterns. Unlike conventional SLMs, in which each pixel must be individually addressed, the proposed SLM scheme generates rich spatial modulation frames under a single uniform optical stimulus, significantly simplifying system complexity and improving stability.

To demonstrate imaging capability, we extract seven representative frames (Fig. 5c and d) from the transient modulation dynamics induced by a single uniform optical pulse to form a modulation matrix with high spatial randomness and orthogonality—key attributes for successful reconstructed images. Using these seven frames, we acquire the corresponding intensity vector I from the SPD and reconstruct binary images representing the target scene using regularization algorithm (see Supporting Information S14). The reconstructed images (Fig. 5f) clearly resolve the target patterns (“T”, “H”, and “Z”), with high fidelity quantified by a MSE that remains close to zero for compression ratios exceeding 66.7% (Fig. 5e). Robustness against system noise is also confirmed, and the reconstruction remains accurate for noise levels up to 0.50% (Supporting Information S16 and S17). This demonstration establishes the hybrid metasurface as a unified platform for ultrafast spectral and spatial modulation, providing a scalable, CMOS-compatible, and optically driven route toward next-generation THz imaging systems.

Methods
Fabrications
Additional fabrication details are provided in Supporting Information S1. The sample substrate consisted of silicon-on-sapphire (SOS), featuring a 460μm-thick sapphire base and an 850nm-thick silicon epitaxial layer. A patterned photoresist layer was first formed on the silicon surface using UV photolithography. A 200nm-thick aluminum film was then deposited on top of the photoresist via thermal evaporation. The sample was subsequently immersed in dimethyl sulfoxide to lift off the photoresist along with the unwanted aluminum, leaving behind the designed aluminum microstructure. Next, UV lithography was performed again to cover and protect both the silicon pattern and aluminum structures with photoresist. The exposed silicon epitaxial layer was removed using inductively coupled plasma (ICP) etching. Finally, the remaining photoresist was cleaned using dimethyl sulfoxide and acetone, resulting in the fabricated hybrid metasurface sample.

Measurements
The OPTP consists of three beam paths, including the detection beam path, the generation beam path and the pump beam path. Firstly, femtosecond pulses with 800nm center wavelength, 100fs pulse width and 1kHz repetition frequency were divided into three beam paths by two beam splitters: the detection beam has a fixed optical range and is directly incident into the ZnTe crystal for detection. The generation beam pumps the ZnTe crystal and generates the THz wave undergoing a delay line to adjust the time of arrival of the THz pulse at the sample. The pump beam also undergoes a delay line to control the time of arrival of the pump pulse at the sample so that varying pump relaxation times of the THz wave are measured. Further details of the OPTP optical setup are provided in Supporting Information S2.

Numerical method
Numerical simulations (transmission spectrum matrix) were carried out using commercially available software (CST Microwave Studio) by a finite-element frequency domain solver with unit cell boundary conditions. The refractive index of silicon was set to n=3.4, while the nondispersive refractive index of sapphire was set to n=3.6. The conductivity of aluminum was defined as 7×106S/m. The eigenmode analysis of the metasurface was carried out in COMSOL Multiphysics.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript makes a strong contribution by reimagining THz-TDS through metasurface-enabled computational spectroscopy and imaging. The integration of ultrafast all-optical modulation with broadband low-correlation metasurfaces is both technically innovative and practically impactful. This work stands out for its conceptual clarity and technical depth, and I identify three key innovations in reshaping the THz-TDS configuration, interesting hybrid metasurface scheme, and novel all-optical ultrafast modulation fashion. Current commercial THz-TDS systems require bulky mechanical delay lines and coherent detection paths for pulse recovery, along with raster-scanning modules for spatial information. This work replace these mechanical elements and coherent detection components with an ultrafast modulated metasurface that encodes spectral and spatial information for computational reconstruction, which represents a substantial simplification of system architecture. The authors present a clever hybrid metasurface design supporting broadband, low-threshold, and low-correlation modulation, based on the theory of bound states in the continuum. Rather than relying on a single narrowband resonance with strong local fields, they employ a band-folding approach to accumulate a dense set of high-Q resonances over a broad bandwidth. This offers a new strategy for tailoring spectral trajectories under a uniform external stimulus. The use of all-optical control provides picosecond-scale modulation speeds without the need for electronic addressing. This not only enhances speed but also highlights the potential for fully integrated ultrafast terahertz platforms. Overall, I consider the work of high quality and recommend it for publication in Nature Communications. Before publication, some points and questions should be addressed by the authors as below:\n1. In terms of systematic performance and resolution, is there a pathway to further improve both spectral and spatial resolution, either through metasurface hardware design or computational reconstruction algorithms? Resolution and reconstruction speed are usually contradictory in conventional setups, but what is the current performance benchmark? I strongly suggest including a table comparing state-of-the-art performance and underlying techniques to clearly situate the contribution.\n2. The manuscript explains that the real parts of the eigen solutions are folded to the gamma point, producing additional resonances. However, it remains unclear how the imaginary parts are controlled. This is essential, as linewidth directly impacts the correlation coefficient.\n3. The hybrid metasurface is fabricated using two photolithography steps with precise alignment on an epitaxial silicon film. Could similar low-correlation modulation be realized on a simpler uniform silicon film deposited by sputtering or related methods?\n4. The ultrafast response relies on fast photocarrier relaxation and large tunability under optical excitation. Why was silicon specifically chosen? Are there material tradeoffs (e.g., carrier lifetime, absorption loss, modulation depth) compared to alternative semiconductors?\n5. In spectral reconstruction, the PSNR improves with more measurement frames, yet the authors only demonstrate reconstruction with 10 frames as the measurement matrix. Why not extend to a larger frame number? What limits the maximum?\n6. In addition, since the photocarrier relaxation is exponential, evenly spaced time slots may not be optimal. How could timeslot selection and frame number be optimized for best reconstruction fidelity?\n7. For reconstruction algorithms, both spectral and imaging reconstruction solve the same inverse problem I=ΦX, yet the algorithms used differ. Elaborate on the rationale, why is a Gaussian basis expansion employed in spectral reconstruction, while imaging uses a different approach? What motivates this distinction?", "reviewer_2": "This manuscript presents an approach to terahertz (THz) spectroscopy and imaging using all-optical hybrid metasurfaces based on quasi-bound states in the continuum (q-BICs). The authors demonstrate a computational spectrometer and a single-pixel imaging system that eliminate the need for mechanical delay lines and coherent detection schemes typically required in conventional THz time-domain spectroscopy (THz-TDS) setups. While the work is technically sound and of interest to specialists in the field, it suffers from several critical limitations that currently restrict its broader impact and general applicability. Given these constraints, I believe the study would be better suited for a more specialized journal. My detailed comments are as follows:\n1. The combination of metallic split-ring resonators (SRRs) and semiconductors (specifically Si and Ge) for investigating ultrafast photoactive dynamics in the THz range constitutes a key aspect of this manuscript. However, a similar approach has already been reported in Laser & Photonics Reviews 2021, 15, 2100244. In that work, the resonant modes—termed the LC mode and dipole mode—exhibit quality factors (Q-factors) comparable to those of the quasi-bound states in the continuum (q-BICs) described here. Given the similarity in both structure and resonant characteristics, it is important to clarify whether these are essentially the same modes.\n2. In the abstract, the authors claim that “achieve dense high-Q resonances across a broad spectral range (0.30-0.55 THz)”. However, a bandwidth of 0.25 THz is very narrow for a practical THz spectrometer, limiting the system's spectral coverage and applicability. Furthermore, the reported Q-factors of the resonant modes are on the order of several tens, which cannot be considered particularly high in the context of modern metasurface-based resonators. For comparison, recent work published in Applied Physics Reviews (1 March 2025; 12(1): 011422) demonstrated q-BICs with experimentally measured Q-factors reaching as high as 371.\n3. In Fig. S6 in SI, the unknown spectral function of random incident light was expanded using Gaussian basis functions. However, Gaussian functions do not form an orthogonal or complete basis set, which raises concerns about the robustness and generality of the reconstruction method. As shown in Fig. 4d, the tested spectral shapes can be easily represented as superpositions of Gaussian functions, which likely contributes to the high quality of the reconstructed results. To better evaluate the performance and versatility of the approach, it is important to test more challenging and non-Gaussian spectral profiles—such as combinations of rectangular or triangular waveforms. Can the reconstruction accuracy be maintained for such cases? Including such tests would strengthen the validity and practical relevance of the proposed method.\n4. As shown in Fig. S7 in SI, even a small noise level of 0.20% can lead to significant distortion in the reconstructed spectrum. This high sensitivity to noise severely limits the reliability and practical applicability of the proposed spectrometer in real-world conditions, where higher noise levels are typically unavoidable. The authors should address this critical issue, for instance, by incorporating noise-robust reconstruction algorithms or regularization techniques, and by providing a more comprehensive analysis of the system’s signal-to-noise performance.\n5. While the authors demonstrate 3×3-pixel imaging of simple characters (\"T\", \"H\", \"Z\"), this resolution is several orders of magnitude lower than what is required for practical imaging applications. Although the authors claim that “the achievable spatial resolution in single-pixel imaging can be improved by designing spatial light modulators (SLMs) with more pixels and optimizing reconstruction algorithms,” the fundamental limitations of their current approach must be addressed. As shown in Fig. 5, each pixel corresponds to a unit array with a distinct resonant frequency. Consequently, even with optimization, the achievable pixel count remains fundamentally limited, making it difficult to reach the resolution levels necessary for real-world imaging scenarios.\n6. The imaging samples used in the demonstration are highly idealized, with transmittance taking only binary values—either 0 or 1—at any given frequency. In practical applications, however, most samples exhibit transmittance values between 0 and 1, and their transmission characteristics often vary significantly across frequencies, potentially including frequency-dependent features such as shifting or broadening of transmission peaks. The current method's performance under such realistic and spectrally complex conditions remains unclear. It is therefore essential to discuss how the proposed imaging approach would handle samples with non-binary, frequency-varying transmittance, and whether the reconstruction algorithm can accurately account for such spectral diversity. Addressing this limitation is crucial for assessing the method’s practical utility.\n7. In Fig. 3g, the spectrum does not fully recover even after 1.9 ns. What is the physical origin of this incomplete recovery?\n8. In Figs. 3b, 3c, 3f, and 3g, there appear to be significant differences in the absolute transmittance magnitudes across the spectra. These variations in signal intensity could influence the accuracy and consistency of the reconstructed results. Why not employ a uniform color scale across all subfigures?\n9. On pages 10 and 11, the authors claimed “the proposed system leverages an active hybrid metasurface as an SLM in conjunction with an SPD, thereby significantly reducing the complexity, size, and cost of the spectrometer”. While it is true that the system eliminates the need for mechanical delay lines and coherent detection schemes used in conventional THz-TDS, it still requires a femtosecond laser-based optical pump system to excite the photoactive metasurface. Such ultrafast laser sources are typically complex, bulky, and expensive—key components that largely determine the overall system footprint and cost. Therefore, replacing one set of sophisticated components (e.g., delay stages and photoconductive antennas) with another (femtosecond laser system) does not necessarily lead to a significant reduction in complexity, size, or cost. This claim should be revised or more carefully qualified to reflect the trade-offs involved, and a more balanced comparison with conventional THz-TDS systems should be provided.\n10. In Fig. 2, the axis labels lack clear indication of values and units, which hinders the quantitative interpretation of the data. Additionally, in Fig. 2a, why is the experimentally measured transmittance consistently lower than the simulated results across the entire frequency range.\n11. On page 3 in SI, the authors present an equation for the frequency-resolved photoconductivity. However, the derivation of this equation is not clearly explained—how was it obtained? Furthermore, the expression appears to account only for the real part of the photoconductivity. What about the imaginary part? The imaginary component is closely related to the dispersive (reactive) response of the material and influences the refractive index of silicon, which in turn affects the resonant frequency, linewidth, and overall behavior of the metasurface modes. Neglecting the imaginary part may lead to an incomplete or inaccurate description of the ultrafast photo-induced dynamics.\n12. On page 3 in SI, how was the equation describing the subsequent relaxation (recombination) dynamics derived? As time increases, the fitted curve approaches a value of approximately –0.032 rather than zero. What is the physical origin of this non-zero offset? Does it arise from a persistent photo-induced change in the material properties, incomplete baseline recovery, or an artifact of the measurement or fitting procedure? Furthermore, in this equation and in Fig. 3, the normalized quantity ΔT/T is used, whereas in Fig. 5, the unnormalized ΔT is presented. What does the parameter T represent in this context? Finally, according to the relaxation equation, the dynamics appear to be independent of the scaling factor S. However, in Fig. 5c, ΔT shows a strong dependence on S. How can this discrepancy be reconciled?", "reviewer_3": "This manuscript presents a compact THz spectroscopic platform that integrates hybrid metasurfaces with computational spectroscopy and single-pixel imaging. The authors introduce a quasi-BIC-based metasurface that achieves ultrafast all-optical modulation, a broad operational bandwidth, and low-correlation encoding. This work directly addresses significant limitations of current terahertz time-domain spectroscopy systems, which depend on bulky mechanical delay lines and raster scanning. The presented platform has the potential to significantly reduce system complexity and cost, thereby accelerating the adoption of THz technologies in practical settings. The use of low-correlated hybrid metasurfaces represents a notable conceptual advance, and the demonstration of computational spectral and imaging reconstruction marks a substantial leap toward practical applications. The experimental demonstrations are robust, providing convincing evidence for high-quality spectral and imaging reconstruction, and the conceptual advance is both original and impactful. Overall, this work reports a compelling scientific advancement, and the manuscript is well-written, timely, and of high significance. I recommend publication in Nature Communications after the authors address the following points:\n1.Major points on benchmarking and comparison with state-of-the-art:\n• A clearer quantitative benchmarking against conventional THz-TDS systems is required in terms of speed, SNR, spectral resolution, and imaging resolution.\n• Provide a direct comparison with competing single-pixel imaging approaches (e.g., compressive sensing using Hadamard or random masks), emphasizing differences in acquisition speed, hardware simplicity, and reconstructed image quality as a function of the number of measurements.\n• The “low spectral correlation” feature of the hybrid metasurface is central to this work. Benchmark this property against existing terahertz modulators, such as graphene, liquid crystal, MEMS, and phase-change metasurfaces, in terms of speed, bandwidth, and encoding quality.\n2.A control experiment is essential to more clearly distinguish the contributions of metasurface design and active material properties. For example, comparing the performance of a silicon film without the quasi-BIC metasurface could elucidate how much of the decorrelation stems from structural engineering rather than the intrinsic material response.\n3.A more thorough discussion of the metasurface design strategy is necessary to clarify how multiple quasi-BICs are managed, how their interference results in low correlation, and why this method is superior to other multiplexing strategies.\n4.The definition of bandwidth should be clearly articulated, such as the -3 dB criterion, and the actual limits of bandwidth should be examined, considering that resonance overlap constrains it. How could bandwidth be extended based on the existing metasurface configuration?\n5.How does encoding efficiency scale with pixel count in single-pixel imaging? A discussion of trade-offs, particularly between modulation depth and pixel density, would be valuable for understanding the limits of achievable spatial resolution.\n6.Additional data on the repeatability, stability, and long-term reliability of dynamic modulation would bolster claims regarding the robustness of hybrid metasurfaces.\n7.How sensitive is the encoding to polarization state and incidence angle? A discussion on robustness would be helpful.\n8.A brief overview of how this platform could be scaled to higher pixel densities and expanded to other frequency ranges would highlight its translational potential. Discuss limitations regarding device area, uniformity, and the potential for extension to other active material platforms.\n9.Some figures (such as figure 1) display low resolution and contrast. It is recommended to provide higher-quality diagrams.\n10.It is recommended to include a coordinate system in Figure 2(a) as the x-axis is referenced several times in the text."}, "author_response": "We would like to thank all the reviewers for their valuable comments, which have helped us improve the manuscript significantly. We have carefully considered each comment and made specific revisions accordingly. Our detailed responses are outlined below.\nPS: Referees' comments are highlighted in 'blue' text, and authors’ responses are shown in 'black' text.\n----------------------------------------------------------------------\nReviewer #1\n----------------------------------------------------------------------\nThis manuscript makes a strong contribution by reimagining THz-TDS through metasurface-enabled computational spectroscopy and imaging. The integration of ultrafast all-optical modulation with broadband low-correlation metasurfaces is both technically innovative and practically impactful. This work stands out for its conceptual clarity and technical depth, and I identify three key innovations in reshaping the THz-TDS configuration, interesting hybrid metasurface scheme, and novel all-optical ultrafast modulation fashion. Current commercial THz-TDS systems require bulky mechanical delay lines and coherent detection paths for pulse recovery, along with raster-scanning modules for spatial information. This work replace these mechanical elements and coherent detection components with an ultrafast modulated metasurface that encodes spectral and spatial information for computational reconstruction, which represents a substantial simplification of system architecture. The authors present a clever hybrid metasurface design supporting broadband, low-threshold, and low-correlation modulation, based on the theory of bound states in the continuum. Rather than relying on a single narrowband resonance with strong local fields, they employ a band-folding approach to accumulate a dense set of high-Q resonances over a broad bandwidth. This offers a new strategy for tailoring spectral trajectories under a uniform external stimulus. The use of all-optical control provides picosecond-scale modulation speeds without the need for electronic addressing. This not only enhances speed but also highlights the potential for fully integrated ultrafast terahertz platforms. Overall, I consider the work of high quality and recommend it for publication in Nature Communications. Before publication, some points and questions should be addressed by the authors as below:\nResponse from Authors:\nWe sincerely thank the reviewer for the careful assessment and positive comments on manuscript, and we appreciate the recognition of our work. Below we provide point-by-point responses to address all comments and concerns.\n1\n1. In terms of systematic performance and resolution, is there a pathway to further improve both spectral and spatial resolution, either through metasurface hardware design or computational reconstruction algorithms? Resolution and reconstruction speed are usually contradictory in conventional setups, but what is the current performance benchmark? I strongly suggest including a table comparing state-of-the-art performance and underlying techniques to clearly situate the contribution.\nResponse from Authors:\nWe sincerely thank the reviewer for this critical and insightful comment. We fully agree that a clear discussion of the potential pathways for improving both spectral and spatial resolution, as well as a benchmarking comparison with state-of-the-art systems, will substantially enhance the completeness and broader relevance of our manuscript.\nFrom a hardware perspective, resolution can be systematically improved through metasurface engineering and system-level optimization. Specifically:\n(i) Enhanced modulation depth: By optimizing the resonator geometry and semiconductor integration to increase the carrier-induced spectral contrast, the modulation depth and resulting signal-to-noise ratio (SNR) can be significantly improved, which directly benefits spectral fidelity and robustness.\n(ii) Low-correlation measurement matrices: Reducing mutual correlation among the encoded spectral or spatial masks can enhance the system’s effective resolution. This can be achieved through dispersion-engineered metasurfaces with denser quasi-BIC modes or through tailored temporal modulation strategies to expand the modulation diversity.\n(iii) Advanced reconstruction algorithms: Incorporating model-based or machine-learning-assisted reconstruction methods can improve both reconstruction speed and accuracy, mitigating the conventional trade-off between resolution and computational cost.\nWe have now included a comprehensive benchmarking summary (Tables R1 and R2 in the revised Supplementary Information S19–S20) comparing our system with representative state-of-the-art THz spectroscopic and imaging techniques in terms of spectral resolution, spatial resolution, acquisition speed, and operational principles.\nThe benchmarking results clearly show that our approach achieves a unique combination of ultrafast, compact, and circuit-free THz spectroscopic and imaging capabilities that bridge the gap between mechanical-delay-free computational methods and metasurface-based on-chip integration.\nWe have added a detailed discussion on potential performance optimization strategies and benchmarking results in the revised manuscript (lines 25~27, page 15 and 3~8, 2 page 16) and Supplementary Information S19–S20.\nAdded S19 Additional comparison of THz spectrometers performance\nTable R1. Performance comparison of common THz spectrometers.\nResolution Peak dynamic Complexity\nScheme Bandwidth (THz) Speed Ref.\n(GHz) range (dB) & cost\nTHz-TDS 0.1-3 ms-s 2-50 60-90 medium S6-9\nCW Photomixing 0.3-3 ms 0.005 - high S10\nFTS 0.1-3 s-min 6 ≈40-50 medium S11\nVNA extension 0.075-1.1 ms-s ≤0.0001-0.01 70-90 high S12\nFilter 0.332-0.377 real-time 0.9-1.0 - high S13\nDispersive Optics 1.3-2.1 s 30-70 - high S14\nTHz-CS 0.3-0.55 ns 30 - low This work\nReference\nS6. Vieweg, N. et al. Terahertz-time domain spectrometer with 90 dB peak dynamic range. Journal of Infrared, Millimeter, and Terahertz Waves 35, 823-832 (2014).\nS7. Nilforoushan, N. et al. Ultra-broadband THz pulses with electric field amplitude exceeding 100 kV/cm at a 200 kHz repetition rate. Optics Express 30, 15556-15565 (2022).\nS8. Bodrov, S.B., Stepanov, A.N. & Bakunov, M.I. Generalized analysis of terahertz generation by tilted-pulse-front excitation in a LiNbO3 prism. Optics Express 27, 2396-2410 (2019).\nS9. Koulouklidis, A.D. et al. Observation of extremely efficient terahertz generation from mid-infrared two-color laser filaments. Nature Communications 11, 292 (2020).\nS10. Hindle, F., Cuisset, A., Bocquet, R. & Mouret, G. Continuous-wave terahertz by photomixing: applications to gas phase pollutant detection and quantification. Comptes Rendus Physique 9, 262-275 (2008).\nS11. Naftaly, M., Dean, P., Miles, R.E., Fletcher, J.R. & Malcoci, A. A simple interferometer for the analysis of terahertz sources and detectors. IEEE Journal on Selected Topics in Quantum Electronics 14, 443-448 (2008).\nS12. Mehdi, I., Siles, J.V., Lee, C. & Schlecht, E. THz Diode Technology: Status, Prospects, and Applications. Proceedings of the IEEE 105, 990-1007 (2017).\nS13. Endo, A. et al. Wideband on-chip terahertz spectrometer based on a superconducting filterbank. Journal of Astronomical Telescopes, Instruments, and Systems 5, 035004 (2019).\nS14. Kanda, N., Konishi, K., Nemoto, N., Midorikawa, K. & Kuwata-Gonokami, M. Real-time 3 broadband terahertz spectroscopic imaging by using a high-sensitivity terahertz camera. Scientific Reports 7, 42540 (2017).\nAdded S20 Additional comparison of THz-SPI systems performance\nTable R2. Performance comparison of common THz-SPI systems.\nImage Compression Complexity &\nModulation mechanism Coding Speed Ref.\nsize (%) cost\nPlanar screen Random - 32×32 0 low S15\nElectrically tuned\nHadamard 1.43 s 8×8 52 high S16\nmetasurface SLM\nDMD Hadamard 10 ms 32×32 0 high S17\nRandom &\nDMD ~ms 128×128 0 high S18\nHadamard\nDMD Walsh-Hadamard ~ms 16×16 0 high S19\nLC metasurface SLM Hadamard ~ms 2×2 50 medium S20\nDMD Hadamard/Walsh - 32×32 0 medium S21\nDSRR array Intrinsic response 2 ns 3×3 22 low This work\nReference\nS15. Chan, W.L. et al. A single-pixel terahertz imaging system based on compressed sensing. Applied Physics Letters 93, 121105 (2008).\nS16. Watts, C.M. et al. Terahertz compressive imaging with metamaterial spatial light modulators. Nature Photonics 8, 605-609 (2014).\nS17. Stantchev, R.I., Yu, X., Blu, T. & Pickwell-MacPherson, E. Real-time terahertz imaging with a single-pixel detector. Nature Communications 11, 2535 (2020).\nS18. Zhao, J., E, Y., Williams, K., Zhang, X.C. & Boyd, R.W. Spatial sampling of terahertz fields with sub-wavelength accuracy via probe-beam encoding. Light: Science & Applications 8, 55 (2019).\nS19. Olivieri, L. et al. Hyperspectral terahertz microscopy via nonlinear ghost imaging. Optica 7, 186–191 (2020).\nS20. Li, W. et al. Dual-color terahertz spatial light modulator for single-pixel imaging. Light: Science & Applications 11, 191 (2022).\nS21. Stantchev, R.I. et al. Noninvasive, near-field terahertz imaging of hidden objects using a single-pixel detector. Science Advances 2, e1600190 (2016).\n2. The manuscript explains that the real parts of the eigen solutions are folded to the 4 gamma point, producing additional resonances. However, it remains unclear how the imaginary parts are controlled. This is essential, as linewidth directly impacts the correlation coefficient.\nResponse from Authors:\nWe thank the reviewer for this insightful comment. We agree that the imaginary parts of the eigenfrequencies, Im[ω], evolving under the band-folding process, indeed play a crucial role in determining the linewidths (radiative losses) and thus the correlation coefficients in the spectral modulation.\nTo clarify this point, we have substantially expanded the discussion in Supplementary Information S3. When the lattice period is doubled along the x-axis, the Brillouin zone (BZ) is folded accordingly, and the modes at the BZ boundary (X point) are mapped to the Γ point. The imaginary components Im[ω] of these folded modes inherit their radiative coupling symmetry from the parent modes. To make this mechanism explicit, we have added Fig. R1c showing the simulated dependence of Q-factors on the wavevector. The results confirm that linewidth tuning can be achieved continuously, enabling flexible control of the resonance decorrelation and modulation depth across the broadband spectral range.\nThese additions clarify how the imaginary parts of the eigenfrequencies are physically engineered, thus directly influencing the observed spectral correlation behavior.\nAdded S3 Band-folding\nAs shown in Fig. R1a, the lattice period of the two DSRRs array is doubled along the x-axis compared with the DSSRs configuration, resulting in a corresponding Brillouin zone (BZ) folding along the Γ–X direction. The outer black square and inner orange rectangle represent the BZs of the DSRRs and two DSRRs lattices, respectively. When the lattice periodicity is doubled, the guided modes located at the BZ edge are folded back to the Γ point (Fig. R1b). The associated radiative Q-factors of these modes follow the same folding behavior (Fig. R1c), while the half band at larger wavevectors in the DSRRs lattice is discarded and reappears as a new band in the folded BZ of the two DSRRs structure. Consequently, the two DSRRs lattice supports three BIC modes that can be coupled to free-space radiation in the interested frequency range. Upon introducing a controlled symmetry perturbation, these BICs evolve into three high-Q q-BIC resonances with distinct mode profiles and enhanced sensitivity, enabling broadband spectral modulation with intrinsically low correlation coefficients.\n5\na b 0.6 c\n108\ny x DSRRs Y(Y′)M Tw ’ o M DSRRs n c y ( T H z) 0 0 . . 4 5 B B I I C C Ⅲ Ⅱ Band folding a ct or 1 1 1 0 0 0 5 6 7 BIC Ⅰ Band folding\nu e Light line Q -f 104\nq\nFr\ne0.3 BIC Ⅰ 103 BIC Ⅱ\nГ(Г′) X’ X 102\n101BIC Ⅲ\nΓ X' X Γ X' X\nWavevector Wavevector\nFig. R1 (a) Illustration of BZ for the band folding operation where X and M points of DSRRs are folded to X′ and M′ points of a two DSRRs, and Г(Г′) point is fixed. (b) Evolution of the real part during band folding (dark open circles: DSRRs; light open circles: two DSRRs). (c) Evolution of the imaginary part during band folding (dark open circles: DSRRs; light open circles: two DSRRs).\n3. The hybrid metasurface is fabricated using two photolithography steps with precise alignment on an epitaxial silicon film. Could similar low-correlation modulation be realized on a simpler uniform silicon film deposited by sputtering or related methods?\nResponse from Authors:\nWe thank the reviewer for this constructive suggestion. To explore the feasibility of achieving similar low-correlation modulation using a simpler process, we experimentally refabricated hybrid metasurfaces in which the silicon layer was deposited by magnetron sputtering instead of epitaxial growth. The process involved:\n(i) patterning the metallic structures on sapphire substrates via standard photolithography and lift-off; (ii) depositing a silicon film by sputtering; and (iii) completing the final lift-off to form the patterned hybrid metasurface (see Fig. R2).\nHowever, optical pump–terahertz probe (OPTP) measurements on these samples revealed no measurable modulation, in contrast to the strong modulation observed in epitaxial Si devices. The key limitation arises from the fundamental material differences between amorphous Si (a-Si) produced by sputtering and epitaxial Si (epi-Si) grown by chemical vapor deposition (CVD):\nStructural quality: a-Si lacks long-range crystallinity and contains a high density of dangling bonds and defects, which act as carrier traps, severely limiting photoconductivity.\nElectronic properties: a-Si exhibits a larger Tauc gap (≈1.6–1.8 eV) and extremely low carrier mobility, leading to weak free-carrier responses even under intense optical excitation.\nEpitaxial Si advantages: epi-Si films are single-crystalline, lattice-matched to the substrate, and possess low defect densities with high carrier mobility and long lifetimes.\nUnder optical pumping near 800 nm (≈1.55 eV), efficient interband excitation produces significant carrier density modulation, yielding deep and broadband changes in 6 terahertz conductivity.\nTherefore, while the overall metasurface geometry can, in principle, be fabricated on sputtered Si, the absence of crystalline order and limited photocarrier mobility in a-Si fundamentally constrain its modulation depth. Realizing comparable low-correlation modulation thus requires a crystalline Si platform (e.g., epitaxial or laser-annealed Si) to ensure sufficient dynamic range and reproducibility in optical-THz coupling.\nFig. R2 Optical image of the fabricated sample and zoomed-in details.\n4. The ultrafast response relies on fast photocarrier relaxation and large tunability under optical excitation. Why was silicon specifically chosen? Are there material tradeoffs (e.g., carrier lifetime, absorption loss, modulation depth) compared to alternative semiconductors?\nResponse from Authors:\nWe thank the reviewer for this thoughtful and important question regarding material selection and tradeoffs. The choice of silicon was guided by a balance between carrier dynamics, modulation depth, fabrication maturity, and system reproducibility, as detailed below.\nSemiconductors suitable for optically driven THz modulation can generally be categorized into two groups:\nIndirect-bandgap materials (e.g., Si, Ge, GaP) typically exhibit longer carrier lifetimes (from several picoseconds to nanoseconds), which can limit modulation speed but provide stable conductivity changes.\nDirect-bandgap materials (e.g., GaAs, ZnTe, InAs) offer faster carrier relaxation (sub-picosecond to few-picosecond timescales), but often at the cost of higher optical absorption, and lower process compatibility with standard CMOS-based fabrication.\nWe ultimately selected crystalline silicon as the modulation platform due to its exceptional material uniformity, low optical loss, mature microfabrication technology, and stable broadband response in THz experiments. Although silicon’s relatively long carrier lifetime is a known limitation for ultrafast modulation, this can be engineered and optimized through techniques such as ion implantation, defect-state introduction, or substitution with low-temperature–grown GaAs (LT-GaAs) to achieve sub- 7 picosecond carrier recovery in future implementations.\nSince this work aims to establish the proof-of-concept of broadband hybrid metasurface modulation, silicon provides a reliable and reproducible basis to validate the physical mechanism without introducing fabrication or material-induced inconsistencies. For clarity, we have included Table R3 summarizing the representative optical and carrier-relaxation parameters of common semiconductors, highlighting the tradeoffs among modulation depth, carrier lifetime, and material loss.\nTable R3. THz response for various semiconductor samples at the excitation wavelengths and fluences.R1-4\nExcitation Fluence Absorbed density Modulation depth Conductivity\nSample\nwavelengths (nm) (μJ/cm2) (photons/cm2) ΔE/E0 (S/cm)\n400 0.95 9.92×1011 0.012 17.0\nSi\n800 1.02 2.75×1012 0.055 0.6\n400 1.27 1.55×1012 0.027 228.3\nGaAs\n800 1.02 2.76×1012 0.173 33.7\nGe 800 0.76 1.78×1012 0.067 47.2\nGaP 400 1.59 1.99×1012 0.007 4.0\nZnTe 400 1.34 1.83×1012 0.016 28.9\n5. In spectral reconstruction, the PSNR improves with more measurement frames, yet the authors only demonstrate reconstruction with 10 frames as the measurement matrix. Why not extend to a larger frame number? What limits the maximum?\nResponse from Authors:\nWe thank the reviewer for this insightful question regarding the number of measurement frames and the observed PSNR saturation. The reconstruction quality in our system is indeed dependent on the number of independent measurements, which determines the rank and conditioning of the measurement matrix. As shown in Fig. R3, the PSNR improves with increasing measurement number in the range of 5–15, but saturates when the frame count exceeds approximately 20. This saturation indicates that additional frames contribute redundant information, as the correlation between adjacent modulation states becomes non-negligible.\nIn practice, the maximum usable frame number is constrained by several factors:\nTemporal modulation window: All frames are captured within the same carrier- 8 relaxation period of the optically driven metasurface (~ns timescale). Extending the frame count effectively oversamples within a limited dynamic window, leading to diminished independence among frames.\nSignal-to-noise ratio (SNR): Increasing frame number without enhancing modulation depth imposes stricter SNR requirements. Excessive frames amplify cumulative noise and can ill-condition the inversion process during compressed sensing reconstruction.\nAlgorithmic efficiency: Beyond the saturation point (≈20 frames), additional measurements provide negligible gains in PSNR while substantially increasing acquisition time and computational overhead.\nTo achieve a balance between reconstruction fidelity and acquisition efficiency, we chose 10 measurement frames, which yield PSNR ≈ 30 dB, conventionally regarded as the threshold for high-quality reconstruction.R5 This trade-off enables robust reconstruction while maintaining real-time measurement throughput.\nWe have clarified this reasoning in the revised manuscript (lines 18-23 of page 12) and emphasized that further improvement may be achieved through adaptive sampling strategies or algorithmic optimization, rather than simply increasing frame count.\nfrequency-shift\nmode-change\nd\nR\n(\nN\nS\nP30\n0\n5 10 15 20 25\nNumber of Measurements\nFig. R3 Reconstruction quality evaluation based on PSNR depending on the number of measuring frames.\n6. In addition, since the photocarrier relaxation is exponential, evenly spaced time slots may not be optimal. How could timeslot selection and frame number be optimized for best reconstruction fidelity?\nResponse from Authors:\nWe sincerely thank the reviewer for this insightful comment. As correctly noted, the photocarrier relaxation in silicon follows an exponential decay, as confirmed by our time-domain fitting in Supplementary Information S5. This means that uniformly spaced temporal sampling does not represent the most information-rich or low-correlated frame selection strategy for spectral reconstruction.\nWe validated this optimization strategy under the same pump fluence (0.48 mJ/cm2), comparing two schemes:\n(i) Time-averaged sampling, in which frames are uniformly spaced in time, and\n(ii) Field-averaged sampling, in which frames are chosen at approximately equal increments of ΔE/E to account for the exponential decay.\n1.0\nTime-averaging\nField-averaging\n0.8\nn\nati o 0.6\nel\norr 0.4\nC\n-\no\nut 0.2\nA\n0.024 THz\n0.0\n0.026 THz\n-0.2\n0.00 0.05 0.10 0.15 0.20 0.25\nΔ¦ (THz)\nFig. R4 Auto-correlation of the transmission matrices obtained by time-averaging (red) and field-averaging (blue).\n7. For reconstruction algorithms, both spectral and imaging reconstruction solve the same inverse problem I=ΦX, yet the algorithms used differ. Elaborate on the rationale, why is a Gaussian basis expansion employed in spectral reconstruction, while imaging uses a different approach? What motivates this distinction?\nResponse from Authors:\nWe sincerely thank the reviewer for this insightful comment. Although both spectral and imaging reconstructions can be expressed as the same linear inverse problem I = 10 Φ·X, the nature of the unknown variable X and the associated prior information fundamentally differ, motivating the use of distinct reconstruction formulations.\n(i) Spectral reconstruction: continuous domain function.\nThe spectral variable represents a continuous function S(ƒ), where the goal is to recover the frequency-dependent spectral amplitude. In this case, the measurement operator Φ maps the continuous spectrum to discrete intensity measurements. Because many THz spectral features (e.g., vibrational or absorption resonances) exhibit approximately Gaussian line shapes governed by the instrument response, we discretize S(ƒ) using a Gaussian basis expansion. This approach offers two key advantages:\na. It provides a physically meaningful, compact representation of spectra, ensuring numerical stability and reduced ill-conditioning.\nb. It aligns the regularization prior (e.g., Tikhonov or nonnegative least-squares) with the natural smoothness and sparsity of spectral signals.\n(ii) Single-pixel imaging: discrete, spatially sparse domain.\nIn contrast, single-pixel imaging aims to reconstruct a discrete spatial distribution rather than a continuous function with a fixed basis shape. Here, the unknown image X does not follow a known analytic form but is typically sparse or piecewise smooth in the spatial domain. Therefore, we employ L1-regularization, which are better suited to enforcing sparsity and preserving sharp edges during image recovery.\nIn summary, the distinction arises because spectral reconstruction leverages a known continuous prior (Gaussian-shaped resonances), whereas imaging reconstruction deals with discrete, structure-agnostic priors (spatial sparsity). These complementary formulations yield more stable and accurate solutions tailored to the statistical and physical nature of each inverse problem.\n11\n----------------------------------------------------------------------\nReviewer #2\n----------------------------------------------------------------------\nThis manuscript presents an approach to terahertz (THz) spectroscopy and imaging using all-optical hybrid metasurfaces based on quasi-bound states in the continuum (q-BICs). The authors demonstrate a computational spectrometer and a single-pixel imaging system that eliminate the need for mechanical delay lines and coherent detection schemes typically required in conventional THz time-domain spectroscopy (THz-TDS) setups. While the work is technically sound and of interest to specialists in the field, it suffers from several critical limitations that currently restrict its broader impact and general applicability. Given these constraints, I believe the study would be better suited for a more specialized journal. My detailed comments are as follows:\nResponse from Authors:\nWe sincerely thank the reviewer for the careful reading and constructive feedback. We have thoroughly revised the manuscript to address the key concern regarding general applicability and broader impact by adding quantitative benchmarks, additional comparisons, and expanded discussions on scope and limitations. These comprehensive revisions substantially strengthen the manuscript’s general significance, demonstrating that the proposed hybrid metasurface platform provides a viable and scalable route toward compact, high-speed, and CMOS-compatible THz spectroscopic and imaging systems. We are grateful for the reviewer’s insightful feedback, which has guided us to significantly improve the clarity, rigor, and broader relevance of our work.\n1. The combination of metallic split-ring resonators (SRRs) and semiconductors (specifically Si and Ge) for investigating ultrafast photoactive dynamics in the THz range constitutes a key aspect of this manuscript. However, a similar approach has already been reported in Laser & Photonics Reviews 2021, 15, 2100244. In that work, the resonant modes—termed the LC mode and dipole mode—exhibit quality factors (Q-factors) comparable to those of the quasi-bound states in the continuum (q-BICs) described here. Given the similarity in both structure and resonant characteristics, it is important to clarify whether these are essentially the same modes.\nResponse from Authors:\nWe sincerely thank the reviewer for raising this insightful point and for drawing our attention to the relevant work reported in Laser & Photonics Reviews (2021, 15, 2100244). That study indeed demonstrates the potential of metal-semiconductor hybrid metamaterials for all-optical modulation in the THz regime, highlighting the community’s growing interest in ultrafast hybrid photonic systems.\nRegarding the cited LPR 2021 study, the resonances described therein, the LC and dipole modes, are intrinsically radiative in nature, exhibiting moderate Q-factors (~5-10 in the THz range). Their coupling gives rise to plasmon-induced transparency (PIT), an analog of electromagnetically induced transparency (EIT) that emerges from destructive interference between two bright modes.R6,7 As shown in Fig. R5, such PIT typically occurs in specific three-level atomic systems, and is limited by their straightforward near-field coupling and shared radiation channels, which leads to simultaneous modulation of both resonances by the integrated semiconductor and thus restricts the formation of multiple, asynchronously addressable resonances with low inter-mode correlation. Consequently, this mechanism is not directly suitable for constructing a broadband, decorrelated spectral basis required for computational THz spectroscopy.\nIn contrast, the double-gap split-ring resonators (DSRRs) used in our design support symmetry-protected modes that are fundamentally non-radiative, forming bound states in the continuum (BICs). These modes are decoupled from free-space radiation channels due to symmetry mismatch, resulting in theoretically infinite Q-factors. By introducing a controlled symmetry-breaking perturbation, we transform the ideal BICs into q-BICs with tunable and experimentally accessible Q-factors in the range of 10-25, chosen as an optimal balance between radiative lifetime and modulation bandwidth. Moreover, the band folding approach applied to our DSRR array enables systematic control over the number and spectral distribution of q-BIC modes, each exhibiting distinct near-field mode profiles and spatial electric field distributions (see Fig. 2b and Fig. S3). These features give rise to spectrally dense, low-correlation modulation patterns, forming the physical basis for our computational spectroscopy and imaging framework.\nFinally, the LPR 2021 reference has been added and appropriately discussed in the revised manuscript. Notably, the polarization-dependent relaxation times (2 ns and 25 ps) reported therein further substantiate our argument that modulation latency can be reduced through judicious selection of semiconductor materials with faster carrier dynamics, reinforcing the future extensibility of our approach.\nThe related discussion has been added to lines 19~25 of page 6 in the revised manuscript.\na b\nwithout coupling\n|1> with coupling\n e\nc nc\nba\n\np A\nbs\nor\n|2>\n|0>\nFrequency\nFig. R5 (a) Schematic of a three-level EIT system. (b) Absorption spectrum showing the EIT effect.\n2. In the abstract, the authors claim that “achieve dense high-Q resonances across a broad spectral range (0.30-0.55 THz)”. However, a bandwidth of 0.25 THz is very narrow for a practical THz spectrometer, limiting the system's spectral coverage and applicability. Furthermore, the reported Q factors of the resonant modes are on the order of several tens, which cannot be considered particularly high in the context of modern metasurface-based resonators. For comparison, recent work published in Applied Physics Reviews (1 March 2025; 12(1): 011422) demonstrated q-BICs with experimentally measured Q-factors reaching as high as 371.\nResponse from Authors:\nWe thank the reviewer for these valuable comments. The term “broad” in our abstract refers to the relative bandwidth with respect to the central frequency, rather than the absolute frequency span. Specifically, our metasurface exhibits a bandwidth of approximately 0.25 THz centered around 0.42 THz, corresponding to a relative bandwidth of ~60%, which is widely regarded as broadband in the THz and photonics communitiesR8. Nevertheless, to avoid possible ambiguity, we have removed the term “broad” in the revised abstract and clarified the context in the main text (line 3 of page 6).\nIn practical applications such as biomacromolecule identification, pharmaceutical characterization, gas spectroscopy, and polymer analysis, the absorption features of target materials are typically confined within narrow windows of a few hundred GHz. Therefore, the current spectral coverage (0.25 THz) already spans the characteristic fingerprints of many representative materials. We have included Table R4 in the revised Supplementary Information S18 to summarize typical absorption bands of THz active materials, showing that our operational bandwidth is well aligned with these spectral 14 regions. Moreover, since the metasurface architecture is scalable by geometric scaling of the resonators, operation at different central frequencies or multiple stacked bands can be readily achieved.\nTable R4. Spectral characteristics of representative materials in the THz band.\nSample Absorption peak (THz) FWHM (THz) Ref.\n0.53 ~0.05 R9\nLactose\n1.37 ~0.1 R9\n0.557 ~0.02 R10\nCHCN\n3\n0.752 ~0.02 R10\nPTFE 6.1 ~0.1 R11\nNylon 6 9.1 ~0.5 R11\nTheophylline Polymorphs 0.96 ~0.1 R12\nArtesunate 1.9 ~0.15 R13\nRegarding the reported Q-factors, we appreciate the reviewer’s comparison with the dielectric q-BIC metasurface in Applied Physics Reviews 12, 011422 (2025), which achieved an experimental Q up to 371. However, we note that the total Q-factor (Q ) total is determined by the interplay between radiative (Q ) and absorptive (Q ) loss rad abs channels following\n1/Q =1/Q +1/Q . R1\ntotal rad abs\nWhile an ideal BIC can in principle support an infinite Q , the experimentally rad measurable Q is ultimately limited by Q , which depends on both material loss and total abs fabrication imperfections. The dielectric resonators in the referenced APR work (based on intrinsic Si) possess negligible absorption, allowing Q to closely approach Q . total rad\nIn contrast, our metasurface employs metallic resonators integrated with semiconductors, which unavoidably introduce finite ohmic losses, thus constraining Q and consequently Q . Nevertheless, our measured Q values (10-23) are relatively abs total high in the reported metallic metasurfaces in the THz regime, consistent with the physical constraints of conductive platforms. Importantly, the moderate Q regime was intentionally chosen to achieve an optimal balance between field enhancement and bandwidth. Extremely high Q-factors, while beneficial for field confinement, would drastically narrow the resonance linewidth, undermining our goal of low-correlation broadband modulation required for computational spectroscopy and imaging (Fig. S8d-f).\nModified text in lines 1~2 of page 6 of the revised manuscript\nThe intrinsic ohmic losses of metallic resonators restrict the increase of Q-factors to a higher value comparable to dielectric metasurfaces.\nModified text in lines 6~8 of page 15 of the revised manuscript\nResonance Q-factors were engineered at appropriate values to balance bandwidth, local field enhancement and spectral diversity, prioritizing low-correlated modulation and computational reconstruction over extreme confinement.\n3. In Fig. S6 in SI, the unknown spectral function of random incident light was expanded using Gaussian basis functions. However, Gaussian functions do not form an orthogonal or complete basis set, which raises concerns about the robustness and generality of the reconstruction method. As shown in Fig. 4d, the tested spectral shapes can be easily represented as superpositions of Gaussian functions, which likely contributes to the high quality of the reconstructed results. To better evaluate the performance and versatility of the approach, it is important to test more challenging and non-Gaussian spectral profiles—such as combinations of rectangular or triangular waveforms. Can the reconstruction accuracy be maintained for such cases? Including such tests would strengthen the validity and practical relevance of the proposed method.\nResponse from Authors:\nWe sincerely thank the reviewer for this insightful suggestion. We originally adopted Gaussian functions as the target spectral profiles because the characteristic absorption fingerprints of biomacromolecules, pharmaceuticals, gases, and polymers typically exhibit Gaussian-like lineshapes under practical conditions. These profiles therefore provide a physically meaningful and computationally stable basis for benchmarking reconstruction accuracy in our target application scenarios.R14\nTo further evaluate the generality and robustness of our reconstruction algorithm, we fully agree with the reviewer’s suggestion and have performed additional tests using non-Gaussian spectra, including rectangular and triangular waveforms with both narrowband and broadband features. As shown in the newly added Fig. R6, the triangular spectra are faithfully recovered with high accuracy (PSNR=30.25 dB for broadband case; 25.23 dB for narrowband case). The rectangular spectra, which contain extremely sharp spectral edges corresponding to high frequencies, show moderate 16 deviations near the edges (PSNR=13.90 dB for broadband case; 17.36 dB for narrowband case), but their principal distributions and bandwidths remain accurately reconstructed. These results confirm that the proposed method maintains good reconstruction fidelity even for non-Gaussian spectral profiles, demonstrating its robustness and versatility beyond Gaussian-type spectra.\na b\n1.2 1.2\nOriginal Original\n1.0 Reconstructed 1.0 Reconstructed\n0.8 0.8\nmission\n0.6\nmission\n0.6\nTransi\n0.4 PSNR=17.36 dB\nTransi\n0.4 PSNR=13.90 dB\n0.2 0.2\n0.0 0.0\n0.30 0.35 0.40 0.45 0.50 0.55 0.30 0.35 0.40 0.45 0.50 0.55\nFrequency (THz) Frequency (THz)\nc 1.2 d 1.2\nOriginal Original\n1.0 Reconstructed 1.0 Reconstructed\n0.8 0.8 PSNR=30.25 dB\nmission\n0.6\nmission\n0.6\nTransi\n0.4 PSNR=25.23 dB\nTransi\n0.4\n0.2 0.2\n0.0 0.0\n0.30 0.35 0.40 0.45 0.50 0.55 0.30 0.35 0.40 0.45 0.50 0.55\nFrequency (THz) Frequency (THz)\nFig. R6 Reconstruction of non-Gaussian spectra. (a, b) Narrowband and broadband rectangular spectra. (c, d) Narrowband and broadband triangular spectra.\nWe have added the above necessary information in the revised manuscript and reproduced it below for the reviewer’s convenience.\nModified text in lines 24~29, page 11 and 1, page 12 of the revised manuscript\nAs a proof of concept, we evaluated the spectral reconstruction performance using a set of THz band-pass filters with different center frequencies and approximately Gaussian passbands. This choice reflects the intended application of our THz-CS, detecting absorption features of biomacromolecules, pharmaceuticals, and gases, whose line shapes are typically well approximated by Gaussian under practical conditions. The attempts to reconstruct rectangular and triangular waveforms can be found in Supplementary Information S13.\n4. As shown in Fig. S7 in SI, even a small noise level of 0.20% can lead to significant distortion in the reconstructed spectrum. This high sensitivity to noise severely limits the reliability and practical applicability of the proposed spectrometer in real-world conditions, where higher noise levels are typically unavoidable. The authors should address this critical issue, for instance, by incorporating noise-robust reconstruction algorithms or regularization techniques, and by providing a more comprehensive analysis of the system’s signal-to-noise performance.\nResponse from Authors:\nWe sincerely thank the reviewer for this important and constructive comment. We fully agree that noise robustness is critical for ensuring reliable spectroscopic performance under practical conditions. In response, we have implemented noise-robust reconstruction strategies, optimized acquisition design, and provided a detailed quantitative analysis of the SNR characteristics of our system, demonstrating reliable spectral reconstruction even at a 1.00% noise level. The revisions are summarized below.\n(i) Noise-robust reconstruction.\nFollowing the reviewer’s suggestion, we incorporated non-negative Tikhonov (L2) regularization and total variation (TV) regularization into the inverse problem, subject to box constraints on the spectral coefficients (0 ≤ α ≤ 1). The regularization parameter was selected using the L-curve criterion. These modifications effectively suppress spurious oscillations and high-amplitude flattening while preserving the correct spectral peak positions. Compared with unconstrained least-squares reconstruction, the regularized solvers significantly increase the reconstruction PSNR and eliminate most false peaks.\n(ii) Acquisition design and frame number\nWe further optimized the modulation matrix and frame selection to improve the conditioning of the inverse problem and increased the number of frames to 12. This enhancement substantially improves the effective HWHM from 0.026 THz to 0.024 THz. Detailed discussions are provided in our response to Comment 6 of Reviewer #1 (page 9~10).\n(iii) Analysis of the system\nTo quantitatively assess the system performance, we added a comprehensive stability analysis of the measurement OPTP system (Fig. R8) by monitoring the time-domain signal of a bare SOS wafer under realistic experimental conditions. We identified three dominant perturbation sources, laser power fluctuations, opto-mechanical imperfections, and environmental disturbances, and analyzed their relative contributions to signal variation.\nThese additions are presented in the revised the Supplementary Information S9, which now includes the detailed system SNR and stability analysis. We believe these enhancements directly address the reviewer’s concern.\nRevised S9: Influence of noise on computational spectral reconstruction\nThe robustness of the computational spectral reconstruction is estimated by adding random noise at a level of 0.15-1.25% to the intensity vector I. Spectra are then 18 reconstructed with our noise-robust solver, which combines non-negative Tikhonov (L2) and total-variation (TV) regularization under box constraints 0 ≤ α ≤ 1; the regularization weight is selected using an L-curve criterion. In addition, we use the optimized frame selection and increase the number of frames to 12, which improves the effective SNR through averaging. Comparisons between the original and reconstructed spectra at varying noise levels are shown in Fig. R7. Low noise (0.25%-0.50%): the reconstructed spectra closely follow the ground truth; both resonance peaks are accurately recovered in amplitude and frequency. Moderate noise (0.75%-1.00%): peak contrast is mildly reduced and slight broadening appears, but peak positions remain correct and no spurious features are observed. High noise (≥1.25%): reconstruction accuracy degrades; while the two characteristic peaks can still be identified, bias in amplitude becomes visible.\na 1.2 b 1.2 c 1.2\nOriginal 0.25% 0.50%\n1.0 1.0 1.0\nTransi\nmission\n0\n0\n0\n0\n.\n.\n.\n.\n2\n4\n6\nmission\n0\n0\n0\n0\n.\n.\n.\n.\n2\n4\n6\nmission\n0\n0\n0\n0\n.\n.\n.\n.\n2\n4\n6\n8\n0.0 0.0 0.0\n0.30 0.35 0.40 0.45 0.50 0.55 0.30 0.35 0.40 0.45 0.50 0.55 0.30 0.35 0.40 0.45 0.50 0.55\nFrequency (THz) Frequency (THz) Frequency (THz)\nd 1.2 e 1.2 f 1.2\n1.0 0.75% 1.0 1.00% 1.0 1.25%\nmission\n0\n0\n.\n.\n6\n8 mission\n0\n0\n.\n.\n6\n8 mission\n0\n0\n.\n.\n6\n0.4\nTransi\n0.4\nTransi\n0.4\n0.2 0.2 0.2\n0.0 0.0 0.0\n0.30 0.35 0.40 0.45 0.50 0.55 0.30 0.35 0.40 0.45 0.50 0.55 0.30 0.35 0.40 0.45 0.50 0.55\nFrequency (THz) Frequency (THz) Frequency (THz)\nFig. R7 Comparison of original and reconstructed spectra with 0.25%, 0.50%, 0.75%, 1.00%, and 1.25% noise levels.\nIn addition, we evaluated the SNR and stability of the OPTP system by measuring the time-domain signals of a bare SOS wafer under realistic conditions (Fig. R8a). We define the time-domain amplitude SNR as SNR=|E |/σ , where |E | is peak amplitude, pk n pk σ is the noise root mean square. The measured SNR is around 34.2 dB. We also n recorded the signal stability for over 10 hours (Fig. R8b). Peak-amplitude variations can be observed with three main sources in the OPTP system:\n(i) Power noise of the femtosecond laser leading to path-wide variations in signal amplitude;\n(ii) Opto-mechanical imperfections from motion in delay stages and inevitable thermal expansion/contraction, that lead to alignment drift;\n(iii) Environmental disturbances, including laboratory vibration and temperature variations.\n19\na b 8\n6\n0.2 6\nd e (\na.\nu.)\n2\n(a.u.)\n-\n-\n0\n0\n0\n0\n.\n.\n.\n.\n2\n1\n0\n1\n56 5 Ti 7 me (ps) 58 d e (\na.\nu.)\n2\n4\nu u\nplit plit\n0\nm 0 m\nA A\n-2\n-2\n-4\n0 20 40 60 0 5 10 15\nTime (ps) Time (ps)\nFig. R8 Time-domain signals of the OPTP system with (a) SOS wafer sample and (b) repeated measurements on the two DSRRs sample within 10 hours.\n5. While the authors demonstrate 3×3-pixel imaging of simple characters (\"T\", \"H\", \"Z\"), this resolution is several orders of magnitude lower than what is required for practical imaging applications. Although the authors claim that “the achievable spatial resolution in single-pixel imaging can be improved by designing spatial light modulators (SLMs) with more pixels and optimizing reconstruction algorithms,” the fundamental limitations of their current approach must be addressed. As shown in Fig. 5, each pixel corresponds to a unit array with a distinct resonant frequency. Consequently, even with optimization, the achievable pixel count remains fundamentally limited, making it difficult to reach the resolution levels necessary for real-world imaging scenarios.\nResponse from Authors:\nWe sincerely thank the reviewer for the thoughtful and constructive comments regarding the spatial resolution of the proposed single-pixel imaging (SPI) scheme. We fully agree that the current 3×3-pixel demonstration primarily serves as a proof of concept to validate the feasibility of ultrafast all-optical modulation and computational reconstruction. In this revision, we have clarified the physical and technological constraints, and proposed a hybrid strategy for extending the method toward real-world applicability.\nAs the reviewer correctly noted, each pixel in our scheme corresponds to a metasurface unit array supporting a distinct q-BIC resonance, which exhibits a unique low-correlated spectral response under optical excitation. In principle, the spatial resolution can be enhanced by reducing the pixel size and increasing the number of unit arrays. In our implementation, the minimal pixel size is determined by the subwavelength period of the metasurface unit cell (94.5 μm along y and 189 μm along x at a central wavelength of 600 μm), which already approaches the far-field diffraction limit. Consequently, the theoretical upper bound of resolution is not set by pixel geometry but by the optical diffraction limit. Operation at higher THz frequencies would further scale down the unit-cell dimension, thereby enhancing the achievable resolution.\nTo move toward practical imaging, we have clarified in the revised manuscript that further improvements can be realized through hardware scaling and algorithmic enhancement by designing metasurfaces with smaller unit cells and more independent modulation channels, and developing optimized compressive-sensing and adaptive reconstruction algorithms to achieve higher compression ratios and better imaging quality.\nMoreover, we propose a hybrid imaging strategy that combines our nanosecond-scale THz-SPI with mechanical raster scanning. In this configuration, each metasurface snapshot captures multiple spatial channels simultaneously, effectively accelerating the acquisition speed by more than an order of magnitude compared with conventional raster scanning under the same field of view. This hybrid approach offers a practical and scalable pathway for extending our ultrafast SPI system to higher-resolution and real-world THz imaging scenarios.\nThe related discussion has been added in lines 18~23, page 15. For the reviewer’s convenience, the revised manuscript is reproduced below.\nAdded text in lines 18~23, page 15 of the revised manuscript\nThe spatial resolution of the proposed scheme is theoretically limited by minimal pixel size of unit cells in the subwavelength range. However, the larger-density elements of modulation matrix would require a higher sensitivity THz detector for accurate intensity readout. For the scenario of large-area object, a hybridized imaging scheme combined with raster scan is necessary, which provides a clear and immediately deployable pathway for real-world imaging scenarios.\n6. The imaging samples used in the demonstration are highly idealized, with transmittance taking only binary values—either 0 or 1—at any given frequency. In practical applications, however, most samples exhibit transmittance values between 0 and 1, and their transmission characteristics often vary significantly across frequencies, potentially including frequency-dependent features such as shifting or broadening of transmission peaks. The current method's performance under such realistic and spectrally complex conditions remains unclear. It is therefore essential to discuss how the proposed imaging approach would handle samples with non-binary, frequency-varying transmittance, and whether the reconstruction algorithm can accurately account 21 for such spectral diversity. Addressing this limitation is crucial for assessing the method’s practical utility.\nResponse from Authors:\nWe sincerely thank the reviewer for this insightful comment. We fully agree that most real-world samples exhibit non-binary, frequency-dependent transmittance characteristics, and it is essential to clarify the applicability of our THz-SPI approach under such conditions.\nWe emphasize that the current demonstration of binary imaging (transmittance values of 0 and 1) is not intended as an idealized simplification but rather as a representative proof-of-concept scenario relevant to practical THz security and nondestructive testing applications. In such contexts, e.g., detecting metallic or highly absorbing objects such as concealed weapons, components, or defects within opaque packages, the transmittance contrast is effectively binary because metals are nearly opaque to THz radiation. Thus, the “0/1” contrast used in our experiment directly reflects realistic use cases in THz imaging.\nFrom a methodological standpoint, our THz-SPI system reconstructs images through compressed sensing at a specific frequency without spectral information, and binary image can be reconstructed from the intensity matrix of single-pixel detector. Nevertheless, the computational spectroscopy module demonstrated in this work (see Figs. 3-4) inherently provides access to the frequency-dependent spectral response.\nTo further demonstrate the capability for dielectric materials that are not completely opaque to THz radiation, we have added experimental and simulated results for dielectric samples that exhibit certain transmission, as shown in Fig. R9 and discussed in Supplementary Information S17. The reconstructed images confirm that our algorithm remains robust for dielectric samples.\nAdded S17 THz-SPI for imaging of dielectric materials\nFor dielectric material whose transmittance is typically non-zero, the THz-SPI can reconstruct the image with good quality as shown in Fig. R9. We examine the cases [0.2, 1], [0.4, 1], and [0.6, 1] to account for dielectric materials with different transmittance, where [A, B] denotes a sample normalized transmittance A and a free-space transmittance B at THz frequencies. Better noise robustness is observed for samples exhibiting larger transmittance contrast.\n22\n[0.2, 1] [0.4, 1] [0.6, 1]\nNoise=0\nNoise=0.10%\nFig. R9 Reconstructed images for transmittance pairs [0.2, 1], [0.4, 1], and [0.6, 1] under Noise = 0 and Noise = 0.10%.\n7. In Fig. 3g, the spectrum does not fully recover even after 1.9 ns. What is the physical origin of this incomplete recovery?\nResponse from Authors:\nWe thank the reviewer for the insightful question. The incomplete spectral recovery at 1.9 ns arises from the intrinsic carrier recombination dynamics of silicon. In lightly doped, low-defect silicon, minority-carrier lifetimes typically span tens of nanoseconds to microseconds. Therefore, partial recovery is expected within the 1.9 ns window accessible in our experiment, which is limited by the maximum delay length of the pump-probe stage. Importantly, although complete recombination is not captured, the key features of the relaxation dynamics including the initial recovery rates and spectral evolution are fully resolved. These features provide sufficient information for the analysis, and the observed partial recovery does not affect the validity of the results or the conclusions of the study.\n8. In Figs. 3b, 3c, 3f, and 3g, there appear to be significant differences in the absolute transmittance magnitudes across the spectra. These variations in signal intensity could influence the accuracy and consistency of the reconstructed results. Why not employ a uniform color scale across all subfigures?\nResponse from Authors:\nWe thank the reviewer for this valuable observation. Small discrepancies in absolute transmittance across the spectra are expected due to practical experimental factors, including non-ideal boundary conditions, finite sample size (edge effects), fabrication imperfections, and minor variations in sample alignment or tilt during measurements.R15 Despite these factors, the overall spectral trends match the 23 simulations closely, ensuring the reliability of the reconstructed results. To improve clarity and facilitate direct comparison, we have updated the color scales in Figs. 3b, 3c, 3f, and 3g in the revised manuscript so that all subfigures use a consistent scale.\nModified Fig. 3 of the revised manuscript\nConductivity (S/m)\na 4800 2400 1200 600 300b c d 1.00\n0.00 Sim. Exp.\nΔE/E\n-\n-\n-\n0\n0\n0\n.\n.\n.\n0\n0\n0\n6\n4\n2 ~27 ps\n@\nf\n0\nre\n.5\nq\n2\nu\nm\nen\nJ\nc\n/c\ny\nm\n-s\n2\nhift\nW/O PUMP q\nq\nq\n-\n-\n-\nB\nB\nB\nI\nI\nIC\nC\nC\nF P Ⅰ\nⅡ\nⅢ\nAuto-Correlation\n0\n0 0\n.\n. .\n2\n5 7\n5\n0 5 Auto-Correlation 0 0 0 0 . . . . 4 4 5 5 6 8 0 2 0.025 Δ0¦. 0(T 0 3H0z . ) 028 0 . T 03 H 5 z\n0.00\n-0.08 0 600 101000 uumm 1200 1800 0.00 0.05 0.10 0.15 0.20\nt (ps) Δ¦ (THz)\nConductivity (S/m)\nΔE/E\ne\n-\n-\n-\n0\n0\n0 0\n.\n.\n. .0\n0\n0\n0 0\n6\n4\n2\n4800\n~30 ps\n2400 1200\n@ 0\nm\n.5\n6\no\n2\n0\nd\nm\n0\ne-\nJ\nc\n/\nh\ncm\nan\n2\nge\n300f\nFano-type Sim.\ng\nW/O PUMP q\nq\nE -\n-\nB\nB\nx I\nI\nF C\nC\np P . Ⅰ\nⅡ\nh\nAuto-Correlation\n0\n0\n0 1\n.\n. . .\n2\n5\n7 0\n5\n0\n5 0 Auto-Correlation\n0\n0 0 0\n.\n. . .\n4\n4 5 5\n6\n8 0 2\n0.020 Δ¦\n0\n0(T.\n. 0H\n0\nz2)5\n0\nz\n.030\nq-BIC Ⅲ 0.00\n-0.08 0 600 100 um 1200 1800 0.00 0.05 0.10 0.15 0.20\nt (ps) Δ¦ (THz)\nFig. R10 Ultrafast measurement matrix with two types of active hybrid metasurfaces with external perturbation. (a, e) Measured relaxation dynamics of active hybrid metasurfaces at a pump fluence of 0.52 mJ/cm2 with microscopic images of the fabricated hybrid metasurface samples.\n(b, f) Simulated and (c, g) experimental dynamic spectral evolution of the two types of hybrid metasurfaces. (d, h) Auto-correlation function in terms of frequency.\n9. On pages 10 and 11, the authors claimed “the proposed system leverages an active hybrid metasurface as an SLM in conjunction with an SPD, thereby significantly reducing the complexity, size, and cost of the spectrometer”. While it is true that the system eliminates the need for mechanical delay lines and coherent detection schemes used in conventional THz-TDS, it still requires a femtosecond laser-based optical pump system to excite the photoactive metasurface. Such ultrafast laser sources are typically complex, bulky, and expensive—key components that largely determine the overall system footprint and cost. Therefore, replacing one set of sophisticated components (e.g., delay stages and photoconductive antennas) with another (femtosecond laser system) does not necessarily lead to a significant reduction in complexity, size, or cost. This claim should be revised or more carefully qualified to reflect the trade-offs involved, and a more balanced comparison with conventional THz-TDS systems should be provided.\nResponse from Authors:\nWe thank the reviewer for the insightful comment. We note that in typical THz-TDS setups, a femtosecond laser is already required as the source for THz generation via photoconductive antennas or optical rectification. In our proposed system, the pump beam for the active metasurface can be easily derived from this same laser source, 24 without introducing additional laser hardware or complexity. Consequently, our approach eliminates the need for mechanical delay lines and coherent-detection electronics, reducing the overall system complexity, size, and alignment requirements compared to conventional THz-TDS.\nTo clarify this point and provide a more balanced perspective, we have revised the manuscript to better describe the trade-offs and advantages of our scheme. The updated text, along with a comparative discussion in Supplementary Information S19, highlights comparison with conventional THz-TDS systems.\nModified text in lines 10~14, page 11 of the revised manuscript\nUnlike conventional systems that rely on mechanical delay lines and coherent detection schemes, the proposed system leverages an active hybrid metasurface as an SLM in conjunction with an SPD, and divides partial femtosecond pulses from transmitter path to drive the SLM, thereby significantly reducing the complexity, size, and cost of the spectrometer.\nAdded S19 Additional comparison of THz spectrometers performance\nTable R5. Performance comparison of common THz spectrometers.\nResolution Peak dynamic Complexity &\nScheme Bandwidth (THz) Speed Ref.\n(GHz) range (dB) cost\nTHz-TDS 0.1-3 ms-s 2-50 60-90 medium S6-9\nCW Photomixing 0.3-3 ms 0.005 - high S10\nFTS 0.1-3 s-min 6 ≈40-50 medium S11\nVNA extension 0.075-1.1 ms-s ≤0.0001-0.01 70-90 high S12\nFilter 0.332-0.377 real-time 0.9-1.0 - high S13\nDispersive Optics 1.3-2.1 s 30-70 - high S14\nTHz-CS 0.3-0.55 ns 30 - low This work\nReference\nS6. Vieweg, N. et al. Terahertz-time domain spectrometer with 90 dB peak dynamic range. Journal of Infrared, Millimeter, and Terahertz Waves 35, 823-832 (2014).\nS7. Nilforoushan, N. et al. Ultra-broadband THz pulses with electric field amplitude exceeding 100 kV/cm at a 200 kHz repetition rate. Optics Express 30, 15556-15565 (2022).\nS8. Bodrov, S.B., Stepanov, A.N. & Bakunov, M.I. Generalized analysis of terahertz generation by tilted-pulse-front excitation in a LiNbO3 prism. Optics Express 27, 2396-2410 (2019).\nS10. Hindle, F., Cuisset, A., Bocquet, R. & Mouret, G. Continuous-wave terahertz by photomixing: applications to gas phase pollutant detection and quantification. Comptes Rendus Physique 9, 262-275 (2008).\nS11. Naftaly, M., Dean, P., Miles, R.E., Fletcher, J.R. & Malcoci, A. A simple interferometer for the analysis of terahertz sources and detectors. IEEE Journal on Selected Topics in Quantum Electronics 14, 443-448 (2008).\nS12. Mehdi, I., Siles, J.V., Lee, C. & Schlecht, E. THz Diode Technology: Status, Prospects, and Applications. Proceedings of the IEEE 105, 990-1007 (2017).\nS13. Endo, A. et al. Wideband on-chip terahertz spectrometer based on a superconducting filterbank. Journal of Astronomical Telescopes, Instruments, and Systems 5, 035004 (2019).\nS14. Kanda, N., Konishi, K., Nemoto, N., Midorikawa, K. & Kuwata-Gonokami, M. Real-time broadband terahertz spectroscopic imaging by using a high-sensitivity terahertz camera. Scientific Reports 7, 42540 (2017).\n10. In Fig. 2, the axis labels lack clear indication of values and units, which hinders the quantitative interpretation of the data. Additionally, in Fig. 2a, why is the experimentally measured transmittance consistently lower than the simulated results across the entire frequency range.\nResponse from Authors:\nWe thank the reviewer for this insightful comment. To improve clarity and enable quantitative interpretation, we have revised Fig. 2 to include explicit axis labels with values and units.\nThe experimentally measured transmittance in Fig. 2a is consistently lower than the simulated results primarily due to the intrinsic material losses of the silicon and sapphire layers, which are not considered in the simulations. Additional contributions arise from practical factors such as finite-size and edge effects, fabrication imperfections, and measurement uncertainties. These experimental realities naturally reduce the transmittance relative to idealized simulations, and the observed resonances remain in good agreement with theoretical predictions. Importantly, these differences do not affect the accuracy or applicability of the system for THz-CS and THz-SPI.\nModified Fig. 2 of the revised manuscript\n26\na x c\ny\nAluminum a frequency-shift mode-change\nw\nl\nd\nSilicon\np 2p\n0.8 0.8\nB B I I C C Ⅲ Ⅱ Sim.\n0.6 BIC Ⅰ 0.6 Exp.\nT0.4 Band folding T0.4 0.8 0.8\nStatic\nDynamic\n0.6 0.6\n0.2 q-BIC Ⅱ 0.2 q-BIC Ⅰq-BIC Ⅱ q-BIC Ⅲ T 0.4 T 0.4\n0.3 0.4 0.5 0.3 0.4 0.5\n¦ (THz) ¦ (THz) 0.2 0.2\n0.3 0 ¦ .4 (THz) 0.5 0.3 0 ¦ .4 (THz) 0.5\nb\nanti-LC parallel-LC dipole-Like\nd\nCorrelation\ncoefficient\nf\nm\nre\no\nq\nd\nu\ne\ne\n-c\nn\nh\nc\na\ny\nn\n-s\ng\nh\ne\nift\n0 4800\nConductivity (S/m)\nFig. R11 Broadband spectral modulation with multiple BICs in hybrid metasurfaces. (a) Underlying physics of multiple BICs interpreted from reciprocal space and realized by period engineering in real space. (b) Mode pattern analysis with surface current and electric field distributions at three q-BICs resonances. The unit cell structure diagram of multi-BICs metasurfaces and key geometrical parameters: p = 94.5 μm, l = 78 μm, g = 4 μm, w = 10 μm, and a = 8 μm. (c) Two types of hybrid metasurfaces by integrating silicon patches at different positions of metallic resonators. External perturbation is introduced by pulse laser with wavelength at 800 nm whose photon energy is higher than bandgap of silicon to inject photocarriers, and thus induces intense modulation of all the three q-BICs. (d) Correlation coefficient evolution at different conductivities of silicon patches.\n11. On page 3 in SI, the authors present an equation for the frequency-resolved photoconductivity. However, the derivation of this equation is not clearly explained—how was it obtained? Furthermore, the expression appears to account only for the real part of the photoconductivity. What about the imaginary part? The imaginary component is closely related to the dispersive (reactive) response of the material and influences the refractive index of silicon, which in turn affects the resonant frequency, linewidth, and overall behavior of the metasurface modes. Neglecting the imaginary part may lead to an incomplete or inaccurate description of the ultrafast photo-induced dynamics.\nResponse from Authors:\nWe thank the reviewer for the insightful comments and agree that a more complete and rigorous explanation of the frequency-resolved photoconductivity is important. In the revised Supplementary Information S5, we now provide a detailed, step-by-step derivation of the retrieval equations (S1–S9) and clarify that the photoconductivity is actually the frequency-resolved complex values, including both the real and imaginary parts, ()()i(). For clarity, we also added an intuitive discussion explicitly linking the imaginary part Δε"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have satisfactorily addressed all of my comments. I am pleased to recommend the manuscript for acceptance.", "reviewer_2": "Although the issues regarding the bandwidth and the number of pixels remain not resolved, my other concerns have been adequately addressed. Therefore, I support it for publication in Nature Communications.", "reviewer_3": "The author has revised the manuscript in response to the reviewers' comments and suggestions. In principle, I concur with the acceptance and subsequent publication of this paper in NC."}, "author_response": ""}]

Decision letter