Near-field probing of the local density of optical states enhanced by bound states in the continuum in nonlocal metasurfaces

Domain: physics
Published: 2025-11-24
Editorial outcome: major_revision

Abstract
Bound states in the continuum (BICs) are optical modes that remain decoupled from free-space radiation. Symmetry-protected BICs in metasurfaces offer powerful means to control light-matter interactions. A key property governing these interactions is the partial local density of optical states (PLDOS), which describes the number of electromagnetic modes available for a photon to occupy at a specific position, frequency and polarization. Here, we employ a terahertz near-field microscope with dual local probes to directly excite and detect quasi-BICs in finite metasurfaces that possess inversion symmetry, corresponding to a symmetry-protected BIC in their infinite counterpart. We observe a strong enhancement of the PLDOS associated with these modes. As the metasurface size increases, the quasi-BIC evolves into a true BIC, with the quality factor diverge while the PLDOS saturates. This establishes an upper limit for enhanced light-matter interaction by BICs. Our findings pave the way for on-chip metasurfaces with maximum light-matter interaction strengths.

Introduction
Bound states in the continuum (BICs) are optical modes that exhibit exceptional spatial localization without any radiation leakage to the far-field, despite their existence within the light cone1,2. In the absence of material losses, BICs have infinite quality factors (Q) while being tightly confined, which provide an ideal platform for enhancing light-matter interactions3,4, and improving the efficiency of photonic devices5–7. To quantify the light-matter interaction in photonic systems, the local density of optical states (LDOS) serves as the most direct indicator. A high LDOS corresponds to regions where optical processes, such as energy exchange or spontaneous emission, are more likely to occur8,9, with the concomitant increase in the decay rate of dipoles located in these regions10,11. Thus, engineering the LDOS by resonant structures through tuning the Q-factor or mode volume (V) significantly influences fundamental processes in these artificial photonic systems12–16. In a resonant system, such as nonlocal metasurfaces, the LDOS scales with Q/V, suggesting that the infinite Q-factor of BICs would lead to a divergent LDOS in metasurfaces and exceptional light-matter interactions17,18. Therefore, the experimental study of BIC-related LDOS is not only fundamentally interesting but also of critical relevance for the development of potential on-chip applications.

The LDOS, ρ(r0, ω), is defined by the number of optical states available for a point dipole to decay. Importantly, the dipole orientation strongly determines the spatial variations in its decay rate, especially in resonant subwavelength systems. Dipoles with different orientations interact with distinct combinations of polarized electromagnetic fields. Therefore, the partial (or projected) local density of optical states (PLDOS) becomes the more relevant quantity to describe light–matter interactions18. For instance, a dipole p aligned along a unit vector ed interacts with modes whose electric fields are polarized in the same direction. At this level of detail, the PLDOS, denoted as ρp(ed, r0, ω), weighs each mode by the electric field projected along the dipole direction ed. The PLDOS capture the directional dependence of the LDOS, enabling the decomposition into contributions from specific polarizations and spatial orientations. Owing to the metasurface formed by an array of resonating scatterers, the electric field generated by a point dipole can be split into two terms: E(r0, ω) = E0(r0, ω) + Es(r0, ω), where E0(r0, ω) is the direct emission from the dipole, and Es(r0, ω) is the scattered field from the resonators19. The PLDOS enhancement, which is equal to the enhanced energy dissipation rate of the dipole, is determined by the ratio of the imaginary components of the scattered and intrinsic electric field amplitudes at the position of the dipole and along the dipole’s orientation8,20:

P/P0 = 1 + Im{p* (r0, ω) · Es(r0, ω)} / Im{p* (r0, ω) · E0(r0, ω)} (1)

Therefore, a direct measurement of the PLDOS enhancement requires determining the complex electric field at the exact position where the dipole source is located, highlighting the significant complexity of these measurements. Due to these limitations, LDOS enhancements are commonly indirectly obtained from measurements of the lifetime or decay rate of excited states21–23. However, this approach is not applicable to BICs, as they do not couple to the far-field continuum. Instead, near-field techniques with specific excitation or detection methods, are required to investigate these dark modes and visualize their spatial field distribution24–26. For example, Dong et al. used electron energy loss spectroscopy to probe the near-field response of optically inaccessible BIC modes in silicon nano-antenna arrays, revealing strong spatial confinement of the associated electromagnetic energy at nanometer resolution26. Hoof et al. used a double THz near-field probe microscope to measure the near field amplitude and phase of BICs and reveal the symmetry protection that suppresses radiation losses25. With the advancement of near-field techniques spanning from the visible to the terahertz (THz) regime, these methods have become widely utilized for controlling and probing electromagnetic waves in engineered photonic and metamaterial platforms27–30. Among them, THz near-field microscopy integrated with dual microprobes offers a powerful approach for investigating diverse photonic phenomena, including spatial coherence and Fourier components of resonant modes25,31. It should be noted that the experimental observation of a pure BIC is not possible, since the required symmetry protection is inevitably broken in real systems by finite size, fabrication imperfections, or probing.

In this work, we present the direct experimental evidence of the enhanced PLDOS associated with quasi-BICs in finite metasurfaces supporting pure BICs in their infinite counterparts. Combined with the theoretical calculation, we find a significant enhancement of PLDOS at the BIC frequency resulting from the collective response of the resonant array, and observe its saturation with increasing array size. The PLDOS enhancement of quasi-BICs in finite systems has thus an upper bound despite the divergence of the Q-factor due to the increase of the mode volume. Our measurements are validated by numerical simulations. Additionally, we reveal the spatial distribution of the PLDOS for different frequencies, highlighting the strong field enhancement and extreme confinement of the this mode to the surface. These results will guide future designs of finite-size metasurfaces with a maximum light-matter interaction strength.

Results
Symmetry-protected BICs are supported by structures that have inversion symmetry2,32. The odd field symmetry of BICs mismatches with the even field of propagating modes in the continuum, leading to a full decoupling of BICs from the continuum33,34. Here, we investigate a metasurface formed by a periodic array of Au rods with a periodicity of p = 300 μm, where one unit cell consists of two identical rods with a separation of d = 80 μm along the x direction, as shown in Fig. 1a. The length and width of each rod are 200 μm and 40 μm, respectively. The metasurface has been fabricated using standard lithography (details on the fabrication can be found in the “Methods” section). Given the designed dimensions, we expect a BIC mode emerging at a frequency of ~0.395 THz as demonstrated in ref. 35. This is supported by a theoretical analysis using a coupled detuned-dipole model, which shows that the imaginary component of the system’s eigenvalue vanishes at this specific frequency.

THz far-field measurements

To investigate the far-field properties of the metasurface and validate the presence of the BIC, we measured the amplitude transmission spectra T(ω, θ) through the sample at various angles θ of incidence with a far-field THz-TDS system (Menlo K15), and normalized them to the transmission of the bare substrate. The THz electric field was oriented parallel to the long axis of the rods (y direction), as depicted in Fig. 1a. At normal incidence, T(ω, 0) does not show any signature of the BIC (Fig. 1b). The broad resonance in the transmission spectrum at 0.495 THz corresponds to the even mode in the array resulting from the in-phase half-wavelength (λ/2) resonances along the long axis of the rods35. A narrow feature appears in T(ω, θ) around 0.365 THz for off-normal incidence (θ = 40°). These findings align with the expectation that incident THz radiation cannot couple to the BIC at normal incidence, rendering it invisible to the far-field incident plane wave. However, the inversion symmetry is broken under oblique illumination and the BIC evolves into a quasi-BIC that can couple to the continuum. The right panel of Fig. 1b shows the transmission at θ = 20°, 30°, 40° normalized by the transmission at normal incidence, highlighting the leaky quasi-BIC resonance. Detailed measurements and simulations of the angle-dependent transmission and the quasi-BIC can be found in Supplementary Information (Fig. S1).

The off-normal excitation breaks the inversion symmetry of the system, causing the BIC to evolve into a quasi-BIC that becomes accessible from the far-field. As the angle of incidence increases, the resonance redshifts and broadens, indicating a reduction in Q-factor. These spectral changes are more clearly observed in the relative transmission spectra, normalized by the transmission at normal incidence, as shown in the right panel of Fig. 1b. The observed redshift and Q-factor reduction are qualitatively consistent with simulation results in Fig. 1c, obtained using the scuff-transmission method in SCUFF-EM (see Methods). In these simulations, the BIC in an infinite array of rod dimers appears at 0.395 THz and exhibits an infinite Q-factor. We note a slight mismatch between the experimental and simulated resonance frequencies. This mismatch arises from idealized assumptions in the simulations with a lossless, infinite, and homogeneous environment, which differ from real experimental conditions (a more detailed discussion is provided in Fig. S1 of Supplementary Information). Despite this discrepancy, the BIC framework remains an effective tool for the design of nonlocal metasurfaces.

THz near-field observation of the quasi-bound states in the continuum

The inversion symmetry of the metasurface prevents BICs from coupling to the continuum of diffraction modes26,36, rendering them invisible in the far-field. Here, we employ two THz microprobes to locally excite and probe the THz near-field response of the metasurface, where the microprobes function as sub-wavelength photoconductive antennas, as demonstrated in ref. 20. The probes locally perturb the metasurface symmetry, transforming the symmetry-protected BIC into a quasi-BIC with a finite radiative linewidth, thereby enabling coupling to and measurement of the otherwise inaccessible BIC. We stress that this local perturbation of the inversion symmetry takes place only at the excitation and detection spots, in contrast to most works in the literature on quasi-BICs, where the symmetry is broken in every unit cell. This makes our systems sufficiently large arrays, where the finite size effects can be neglected, the closest possible experimental realization and characterization of a pure BIC. Due to the much smaller size of the photoconductive antennas compared to the wavelength of THz radiation (~λ/38 at the BIC frequency), the microprobe acting as a THz emitter can be approximated as a radiating point dipole. The other photoconductive antenna serving as the detector is placed at a fixed distance of 10 μm from the emitter, as shown in Fig. 2a. The details on alignment of microprobes can be found in Supplementary Information (Figs. S2–S4). This close proximity ensures that the emission and detection occur at the “same position” within a sub-wavelength distance as shown in (i) and (ii) of Fig. 2a. In the experiment, both THz antennas are oriented along the y-direction or the x-direction (along the long or short axis of the rods), fixing the orientation of the dipole moment of the THz emitter and the detected field component of the near-field. Therefore, the detected signal is proportional to the PLDOS, Gyy(r0, r0, ω) or Gxx(r0, r0, ω) (see next section), which represents the components of the dyadic Green function along the y- or x-directions at the position of the emitter, respectively. Based on this unique setup, we have previously investigated the PLDOS enhancement above gold, InSb, and quartz interfaces, achieving excellent agreement with theoretical calculations for a point dipole as a function of height and demonstrating that the microprobes do not alter significantly the PLDOS20. In what follows, we simplify the notation by not expressing the spatial and frequency dependencies of the Green function.

Both microprobes are positioned at a height of 10 μm above the surface for spatial mapping of the PLDOS. The time-dependent THz electric field of the locally excited metasurface is spatially mapped on an x-y grid with a step size of 4 μm. While the line spread function of the near-field probes has been measured to be 20 μm (see Supplementary Information, Fig. S5), we chose a smaller step size for a detailed mapping. Figure 2b displays the distribution of the measured THz electric field normalized to the maximum across a unit cell, obtained from the peak value of time-dependent THz amplitude. As shown in this figure, the normalized THz electric field along the y-direction is enhanced at the top and bottom edges of the rods.

To analyze the near-field measurements in more detail, we extract the THz transients at each position. Fig. 2c displays the time-dependent THz electric-field amplitudes at the positions marked in Fig. 2b, which are located at the top edge of the rod, at the corner, and at the edge of the unit cell, respectively. The THz transient measured at the edge of the rod exhibits a pronounced oscillatory behavior with a long decay time of over 35 ps and an increased peak intensity. These oscillations are also present, although less pronounced, in the measurement at the corner of the unit cell. Fourier transformation of the THz transients provides the frequency-dependent intensity of the THz field |E(r0, ω)|2, plotted in Fig. 2d for the three locations. The near-field measurements provide direct evidence of a greatly enhanced electric field at the edge of the rod, observed as the peak in the THz intensity centered at 0.395 THz. This peak, corresponding to the oscillations in the THz transients, originates from the field enhancement owing to the quasi-BIC mode. The peak has a finite linewidth, which results from material losses in the Au rods and the limited spectral resolution caused by the short time window of the THz transients measurements37,38. Additionally, a high but relatively lower field enhancement is observed at the corner of the unit cell shown in Fig. 2d, while it is nearly absent at the edge where the quasi-BIC’s electric field is minimal. The enhancement of the field at the corner of the unit cell illustrates the interaction between neighboring unit cells in the array originating from the THz wave scattering by the Au rods, and highlights the importance of the extended metasurface for the formation of the BIC.

Enhancement of the partial local density of optical states

The dyadic Green’s function, G, is a second-rank tensor represented by a 3×3 matrix for the x-, y-, and z-field components. The trace of the imaginary component G at the position of the source accounts for the total LDOS corresponding to the three possible orthogonal dipole orientations. The PLDOS is used to describe the LDOS of a dipole oriented along a principal axis18,20. The THz microprobes used in the setup generate and detect linearly polarized broadband THz radiation39. Therefore, the setup is particularly suited for measuring the PLDOS.

To gain insight into the PLDOS originating from a quasi-BIC in a metasurface of finite size, we calculated the frequency-dependent PLDOS enhancement with respect to vacuum, Im(Gyy(r0, r0, ω)) / Im(G0(r0, r0, ω)), for a point dipole oriented along the y-direction. The calculations were performed using the package SCUFF-EM, as detailed in the Methods section. Figure 3a, b demonstrates an enhancement near the top and bottom edges of the Au rods for both a single dimer and an array of 9×9 dimers, at their respective resonance frequencies. Moreover, this enhancement increases significantly with the number of unit cells in the array. Unlike local metasurfaces, where each unit cell or meta-atom responds independently to the incident field, nonlocal metasurfaces support optical modes that arise from coherent interactions among many identical neighboring units. This nonlocality is evident in the behavior of the quasi-BIC mode in the metasurface system studied here. Closer inspection of the frequency-dependent PLDOS enhancement at the position of the rod’s edge (see Fig. 3d) reveals a number of peaks around the frequency of the BIC mode for the infinite array, which increase in number and become sharper as the array size increases (Further details are presented in Fig. S6 of Supplementary Information.). These peaks are collective hybridized modes, formed by coherent coupling of localized resonances across the unit cell array. They manifest as standing wave patterns extending over the entire structure. The inset in panel (ii) of Fig. 3d shows a magnified view of the five distinct modes for the 9×9 array. A spatial map of the first three modes is shown in Fig. 3e. (The PLDOS distribution of other modes is presented in Fig. S7 of Supplementary Information.) These modes experience radiative losses due to the finite nature of the array. For a perfect, lossless and infinite system, the collective modes would resonate with infinite Q-factor, defining BICs with zero linewidth40,41.

The THz near-field microscope allows us to experimentally determine the PLDOS enhancement along the same principal y-axis of the metasurface containing 84×84 unit cells. A bare quartz substrate, identical to the sample substrate, was used as the reference to normalize the measurements and isolate the optical response of the metasurface. As illustrated in Fig. 3c, we observe a significant PLDOS enhancement at the frequency of the quasi-BIC in the measurements, in agreement with the calculated results. Fig. 3d, panel (i), shows the measured frequency-dependent PLDOS enhancement at the edge of a central rod in the array. Rather than a series of sharp peaks compared to the results from calculations in panels (ii)–(v), we observe a broad peak in the measured results with a very large enhancement. This difference is due to the limited spectral resolution set by the range of the time-delay stage, which is used to acquire the THz transients. In addition, the differences between the measured and calculated values for each sample could be partly attributed to intrinsic material losses within the samples and refractive index mismatches at the interface. Peaks slightly broader than those in the simulations can be observed in the measurements of the PLDOS enhancement of arrays with sizes of 3×3 and 5×5 dimers, as illustrated in Fig. 3d, panels (iv) and (v). In addition, a very broad resonance with low PLDOS enhancement can be observed in the measurements and calculations of the single dimer (Fig. 3d, panels (vi)). This dependence of the peak width on the array size illustrates the nonlocal resonance of the metasurface due to the collective response of the dimers16,42.

An intriguing question is why the calculated and the experimental PLDOS do not show a diverging behavior at the BIC frequency. A symmetry-protected BIC mode has an infinite lifetime and exists only in ideal systems that extend infinitely in at least one direction. In practice, the BIC arises from the collective resonances of constituent scatterers-in this case, individual dimers26,35,43,44. However, the finite size of the arrays used in both calculations and experiments leads to radiation leakage at the boundaries, resulting in a finite lifetime; thus, the ideal BIC manifests as a quasi-BIC whose lifetime depends on the array size41. Still, the question remains about the relation between the BIC lifetime and the PLDOS. We have calculated the Q-factor, defined as the reciprocal of the mode’s lifetime, as a function of the array size. The Q-factor was extracted from the first peak (fundamental resonance) of the frequency-dependent PLDOS enhancement, which is described in Fig. S8. Figure 3f, g show the calculated Q-factor and the PLDOS enhancement as a function of the total number of dimers in the array, respectively, indicating that while the Q-factor of the quasi-BIC diverges when the size of the metasurface is increased (as theoretically pointed out for finite sphere arrays in ref. 41), the PLDOS enhancement remains finite. The measured values can be found in Fig. S9. Intuitively, the PLDOS depends on the Q-factor and the mode volume V through the relation ~Q/V. While the Q-factor increases with the size of the metasurface, the mode volume V also grows, which explains the saturation (or very slow growth) of the PLDOS enhancement for large metasurfaces despite the diverging Q-factor.

Spatial distribution of partial local density of optical states

Given that the PLDOS enhancement at the frequency of the quasi-BIC mode is supported by collective resonances in the dimers coupled across the array, it is crucial to investigate the spatial distribution of the PLDOS in the metasurface. This distribution defines the relevant mode volume of the quasi-BIC and the extent to which light-matter interactions can be modified by the metasurface45,46. The calculated PLDOS enhancements are normalized to the corresponding values in vacuum at each frequency, whereas the experimental results are normalized to the PLDOS measured on the bare quartz substrate at the same frequency. From the calculations of the 9×9 dimer array, we obtained the spatial distribution of the PLDOS enhancement in the central unit cell and in the xy-plane at a height of 10 μm above the metasurface, as plotted in Fig. 4a–c for 0.35 THz, 0.395 THz, and 0.45 THz, respectively. The enhancement of the PLDOS is strongly localized at the edges of the rods and at the frequency of the BIC for the infinite array (see Fig. 1c), 0.395 THz, rapidly decreasing beyond these areas and frequency. The greatest enhancement occurs in an area of 40×20 μm2, which is almost three orders of magnitudes smaller than λ2. This result illustrates the extremely high degree of spatial confinement of the electric field with a small volume resulting from the origin of the BIC on the λ/2 resonances in the rods25. This strong confinement may even be surprising, as the BIC mode is supported by collective resonances in the dimers coupled across the infinite array. For frequencies slightly away from the BIC frequencies, the small PLDOS enhancement also rapidly decreases with the spatial extent.

Using the near-field microscope, we experimentally determined the spatial extent of the PLDOS by scanning the THz microprobes with a resolution of 4 μm across the xy-plane above the metasurface. We also observe a PLDOS enhancement that is highly confined to the edges of the rods, as shown in Fig. 4d–f. The measured PLDOS enhancement appears more pronounced at the upper ends of the rods than at the lower ones, which is caused by a small intentional misalignment between the THz emitter and receiver, as shown in Fig. 2a. This misalignment prevents cross-talk between emitter and receiver that would be caused by specular reflection of the laser radiation used to gate the photoconductive switches. In addition, the phase information between the two bars within a dimer has been demonstrated and discussed in Fig. S10 of the Supplementary Information. For x-polarized excitation and detection, which determines Pxx/P0, both calculations and experiments indicate that the PLDOS enhancement occurs along the long sides of the rods (cf. Fig. S11 of Supplementary Information).

Out-of-plane confinement of the enhanced partial local density of optical states

In addition to the in-plane spatial distribution of the PLDOS discussed so far, the high confinement of the quasi-BIC should also extend to the out-of-plane direction, confirmed in Fig. S12 of Supplementary Information, which reflects the evanescent decay of the electric field enhancement toward the continuum47. Experimentally, we determined this PLDOS confinement to the surface by positioning the microprobes at (90 μm, 250 μm) for y-polarization and (180 μm, 210 μm) for x-polarization, corresponding to the position in Fig. 2b, and gradually moving them in the vertical (z) direction away from the metasurface with a step size of 1 μm.

Figure 5a shows the THz transients for y- and x-polarized THz excitation and detection, respectively, in (i) and (ii). At 5 μm above the surface, the THz transients exhibit a strong under-damped oscillation, corresponding to the long lifetime of the quasi-BIC. As the distance from the surface increases, this oscillation rapidly weakens and becomes damped. At 35 μm above the surface, the oscillation is barely visible, illustrating the strong field confinement of the quasi-BIC to the surface47.

To clearly show the decay of the quasi-BIC from the surface, the amplitude spectra are obtained by performing a Fourier transform of the THz transient signal, shown in Fig. 5a, after subtracting the contribution from direct dipole interaction. A clear peak at 0.395 THz, which corresponds to the quasi-BIC, is observed in the frequency spectra for both y- and x-polarized THz excitation and detection, respectively, shown in Fig. 5b, and the intensity of the peak rapidly decreases with increasing probe height above the surface.

The tight confinement of the quasi-BIC to the surface is even more clearly visible in the frequency-dependent PLDOS enhancement, which is experimentally obtained by Fourier transforming the THz transients measured on the metasurface and dividing them by the THz transients measured on a quartz substrate, as shown in Fig. 5c, d for both polarizations. The PLDOS enhancement reaches a maximum at the BIC frequency of 0.395 THz for both polarizations and is confined within a height of 20 μm from the surface. The equivalent calculated PLDOS enhancements, Pyy/P0 and Pxx/P0, are shown in Fig. 5e, f and confirm the experimental results.

The spatial extent of the quasi-BIC mode in the xy-plane can be further quantified by examining the normalized PLDOS at the frequency of the quasi-BIC mode as a function of the relative height, h/λ0, as displayed in Fig. 5g. The enhancement of the PLDOS in such a metasurface occurs only within a distance of λ0/20, demonstrating strong out-of-plane confinement of the quasi-BIC. This high spatial confinement, evidenced by the PLDOS enhancements (shown in Figs. 4 and 5), provides strong evidence of the potential of symmetry-protected BICs in metasurfaces to enhance light-matter interaction. This can be leveraged for highly integrated on-chip applications such as sensing, nonlinear optics, and lasing.

Methods
Sample fabrication
The gold metasurfaces were fabricated on a 2.5×2.5 cm2 quartz substrate using a standard lithography process. The AZ5214 UV photoresist was spin-coated on the substrate with a thickness of 2 μm. The exposure was done with a hard mask by UV lithography (MA6 Mask Aligner). Subsequently, the exposed sample was developed in an AZ 726 MIF solution for 1 min. The metal film was prepared by the electron beam evaporation of a 5-nm-thick titanium film as the adhesion layer and an 80-nm-thick gold film. The lift-off process of the samples was done with an acetone solution and the samples were cleaned with ultrapure water.

PLDOS calculations
Numerical calculations of the PLDOS were carried out with the Surface Current/Field Formulation of Electromagnetism (SCUFF-EM), which is a free, open-source software implementation of the boundary-element method (or the “method of moments") of electromagnetic scattering48. In particular, we make use of the tool SCUFF-LDOS for computing the electromagnetic LDOS at points inside or outside compact or extended material bodies. Actually, it computes the scattering part of the dyadic Green’s functions of the geometry in question. For nonperiodic geometries, SCUFF-LDOS does six scattering calculations for each (x, ω) point – specifically, in which the incident field is the field of an electric or magnetic point source oriented in each of the three Cartesian directions. The same geometrical parameters as in the experimental measurements were used, but gold rods were considered as perfectly planar rectangles embedded in a uniform medium with n = 1.55. This refractive index is the average of those of air and the supporting quartz substrate, n = 2.1. This approximation is known to yield consistent results except for very close to the interface49.

THz near-field emission and detection setup
The THz near-field setup consists of two commercially available THz near-field microprobes (TeraSpike TD-800-X-HR-WT, Protemics GmbH). The tip of the microprobes is made of low-temperature-grown GaAs (LT-GaAs), with two microstructured electrodes that define a photo-conductive antenna. A short laser pulse with a duration of 100 fs and a central wavelength of 780 nm is used to generate and probe the THz amplitude with a delay stage as a function of time, i.e., the THz transient pulse, by exciting the electrically biased emitter and the unbiased detector, respectively, as shown in Fig. 2a. Its working principle is the same as the standard far-field THz-TDS, which is a phase-sensitive technique that retrieves the transient THz electric field E(t) as a function of time t50. By applying a Fourier transform to E(t), the complex THz electric field E(ω) can be obtained across a broad frequency range. The frequency resolution of the terahertz near-field time-domain spectroscopy (THz-TDS) system is approximately 20 GHz, determined by the 50 ps time window from the delay stage. The size of the THz source, defined by the emitter microprobe, is ~20 μm, which is much smaller than the period of the metasurface (300 μm) and the wavelength at the frequency of the quasi-BIC (λ0 = 760 μm). Therefore, we approximate the source as a point dipole.

The discriminating advantage of polarization-sensitive microprobes over a standard far-field setup lies in their ability to measure the field component E(t) along the orientation of the detector microprobe at the position r0 where it was emitted (within a subwavelength distance given by the separation between the emitter and detector microprobes), thus allowing the determination of the fully vectorial complex field E(r0, ω)39. Such a near-field THz-TDS system is ideally suited to experimentally determine the PLDOS, which corresponds to the principal components of Im(G(r0, r0, ω)), defined by the orientation of the dipole moment of the emitter microprobe20.

The PLDOS of the metasurface is spatio-temporally resolved at every point via a raster scan (step size: 4 μm) with sub-THz-cycle temporal resolution. Fast Fourier transforms of the THz transients at each position provide the complex spectral information with the real and imaginary components of the THz field. Owing to the limited scan range of the time-delay stage, zero padding is applied to the time-domain data prior to the Fourier transformation to enhance the frequency spectral resolution in the frequency domain. The measured PLDOS enhancements Pyy/P0 and Pxx/P0 are retrieved from the imaginary field component of the measured THz electric near field on the array, Im(E(r0, ω)), under y- or x-polarization of the near-field emitter and detector, normalized by the measurements on the underlying quartz substrate, Im(E0(r0, ω)).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript presents a study on the near-field terahertz probing of metasurfaces supporting bound states in the continuum (BICs), focusing on the enhancement of the partial local density of optical states (PLDOS). The experimental work is original and commendable, and the results are significant. The data and methodology are of high quality, incorporating simulations and experiments, which agree well – the agreement is qualitatively good and the reported effect is convincingly observed - already quite a challenge in these sorts of systems. The data is well presented, and the language is clear.\nHowever, several aspects related to data interpretation, presentation, and contextualization require attention to further enhance the manuscript's clarity and completeness.\n1. The manuscript states, \"we have developed a near-field THz time-domain spectrometer…\" suggesting the development of a new setup. However, unless I am mistaken, a rather similar system has been previously reported by the authors (e.g., ACS Photonics, 2022, 9, 5, 1547–1554). It would be more accurate to acknowledge that this established setup is now being utilized (or modified) to probe BICs in metasurfaces.\n2. The concept of PLDOS is central to this study. It first appears in the abstract without explanation (whereas the LDOS is introduced well). Given its nuanced nature, a clear definition accompanied by relevant literature would greatly benefit readers. Unless I am mistaken, the current manuscript lacks any references on PLDOS, which should be addressed. For example, I note: Ly, Ousmane, et al. \"Partial local density of states from scanning gate microscopy.\" Physical Review B, vol. 73, no. 23, 2006, p. 235305. DOI: http://doi.org/10.1103/PhysRevB.73.235305.\n3. The manuscript mentions “a fixed 10 μm distance between the emitter and detector, as shown in Figure 2(a)”. However, this distance is not clearly depicted in the figure. In my experience, determining and maintaining such spacing in near-field experiments can be challenging. On a related note, I see that experimental data shown in Figure 5c appears to originate at 0μm. Could the authors elaborate on how this distance was measured or verified, including any associated uncertainties?\n4. Simulation results are placed in the supplementary material and seem to exhibit narrower features compared to experimental data. This discrepancy should be elaborated upon. Where possible, I would recommend to:\n- Integrate key simulation results into the main figures – particularly Figure 1b, where the most important feature takes up only a small portion of the figure.\n- Overlay experimental and simulation data on the same axes for direct comparison (e.g., in Figure 3d), as independently scaled axes and the separation into insets may obscure important differences.\n- Apply color scales that do not obscure data, for example, Figure 4 appears to contain saturated / truncated colourmap data, and seems to have been interpolated instead of showing the actual data in each measurement pixel.\n- Each subfigure in Figure S1 b,c, have different frequency axes.\n- A version of Fig. S1c using experimental and the same frequency axes would be beneficial (the data is already there.) Alternatively, a version of Fig. S1b with simulation data might help comparisons (but is less effective at capturing the emergence of coupling to the BIC mode).\n5. The study employs a mapping step size of 4 μm, while to my knowledge the photoconductive antennas typically have an effective resolution of 10 μm according to the specifications. Could this lead to oversampling and feature blurring? Do the authors have an estimate of the system's point spread functon(PSF), or can they comment on the effective resolution?\n6. Figure 3 illustrates the trend of increasing BIC behavior with metasurface size. However, experimental data for N = 7 and N = 9 are absent. Have measurements been conducted for these configurations? Additionally, could signal processing techniques, such as zero-padding, enhance the visibility of spectral features, perhaps in smaller metasurfaces where resonances are wider (e.g., N = 5)?\n7. Are there simulations corresponding to Q and Pyy/Po in Fig. 3(f,g)?\n8. I suggest making a small note clarifying that quartz is used as a reference due to it being the substrate material.\nFurther comments on potential related literature:\nThis work sits within a broader research field using near-field antennas to probe photonic structures (wherein the present paper considers the case where two antennas are very close together). The authors are therefore encouraged to consider whether any of the following works could be relevant and might merit citation:\n- Gansch et al., Light: Sci. Appl., 2016. DOI: https://doi.org/10.1038/lsa.2016.147\n- Yermakov et al., Phys. Rev. B, 2018. DOI: https://doi.org/10.1103/PhysRevB.98.195404\n- Bulu et al., New J. Phys., 2005. DOI: https://doi.org/10.1088/1367-2630/7/1/223\n- Li et al., Phys. Rev. Applied, 2025. DOI: https://doi.org/10.1103/PhysRevApplied.23.014084\n- Lezhennikova et al., Comm. Physics 2023, https://doi.org/10.1038/s42005-023-01347-1\n- Lezhennikova et al., ACS Photonics 2023, https://doi.org/10.1021/acsphotonics.3c01312\nAt this point I would like to declare that I am an author of the last paper of the above list (Lezhennikova et al., ACS Photonics 2023). It presents a related approach using the same combination of antennas (albeit at different locations on the opposite side of a photonic crystal sample), and it may offer complementary insights, particularly regarding a spatial Fourier-space analysis and the observation of enhancement / suppression of a near terahertz field emitter as a result of the available photonic states. Of course I am not requesting that this work be citated in the revised manuscript, but merely making the authors aware of it in case they happen to consider it relevant or useful.\n9. I think it would be broaden the appeal of this work to speculate more, in the conclusions, about the implications of this work at other frequency ranges, or for photonic- and quantum- technologies in general, e.g., in the context of quantum or single photon emitters (if appropriate), and beyond, at the authors’ discretion.\nMinor Suggestions:\n- The measurement conditions (“4 μm resolution, 10 μm height”) are repeated more than once.\n- Clarify the method for Q-factor extraction related to the phrase “The Q-factor was extracted from the first peak (fundamental resonance) of the frequency-dependent PLDOS enhancement”\n- “We have calculated and the Q-factor…” – is this a typo?\n- Please specify the model of the Protemics GmbH antenna used, as there are several and they vary in characteristics.\nIn conclusion, this study presents impressive results in a technically challenging field. The observations are significant, and the technique holds promise. My suggestions aim to improve clarity, transparency, and context. Congratulations to the authors on the experimental achievement.", "reviewer_2": "In the manuscript “Near-field Probing of the Local Density of Optical States Enhanced by Bound States in the Continuum in Nonlocal Metasurfaces” by Jie Ji et al, the authors demonstrate the near-field mapping of BICs in metasurface and vertify the associated PLDOS enhancement. Unfortunately, I find that this manuscript demonstrates some overlaps with the prior works (ter Huurne, Stan, et al. \"Direct Observation of Lateral Field Confinement in Symmetry‐Protected THz Bound States in the Continuum.\" Advanced Optical Material 11.6 (2023): 2202403; van Hoof, Niels JJ, et al. \"Unveiling the symmetry protection of bound states in the continuum with terahertz near-field imaging.\" ACS photonics 8.10 (2021): 3010-3016.). As a result, the results in this manuscript are not particularly innovative and impressive. Regrettably, I cannot recommend this manuscript for publication in Nature Communications. Additional comments are as follows:\n1. In the work, the authors mentioned, “To investigate the far-field properties of the metasurface and validate the presence of the BIC, we measured the amplitude transmission spectra T(ω, θ) through the sample at various angles of incidence θ with a far-field THz-TDS system.” They demonstrate the existence of BIC through angle-resolved far-field experiments. However, the limited comparison of experimental results under only two angular incidence is not sufficient. More importantly, rigorous validation criteria for BIC require explicit verification of the presence of polarization topological charges.\n2. The authors assert direct observation of a BIC mode through near-field measurements.\nWhile the microscopic probe as an excitation source in Figure 2(a) indicates a tilted orientation, which likely excites a quasi-BIC rather than a genuine symmetry-protected BIC. If the tilt angle of microprobe emitter is systematically varied, would the measured radiative leakage (e.g., resonance linewidth and Q-factor) show a significant variation—a hallmark of quasi-BIC behavior?\n3. In Figure 3(d), significant discrepancies are observed between the experimental and simulated curves, including variations in linewidths and resonance frequency shifts. Particularly for N = 5, anomalous multi-peaks are in the measured results. The authors should explain these differences in detail. Although the authors attribute peak broadening to the system’s frequency resolution limitations, they have not provided critical parameters about the delay window and the spectral resolution of their system. In addition, the experimental spectra in Figure 3(d) exhibit insufficient sampling points (the inflections in the experimental curve are overly pronounced). Since the signal fluctuation is also influenced by noise during measuring, this effect becomes particularly pronounced when the number of sampling points is limited.\n4. Regarding the Q factor variations in Figure 3(f), the authors claim that the Q factor is defined as the reciprocal of the mode’s lifetime. They should provide the detailed calculated method of the Q factor and display the experimental trend of Q-factor evolution with the number of dimers (N).\n5. In Figure 5(a) and (b), the calculated Pyy/P0 exhibits intensity distribution below the BIC frequency, whereas the experimental results show prominent intensity above the BIC frequency. The authors should explain this discrepancy. Additionally, the authors should provide frequency-domain spectra at multiple height levels （the results of four height in Figure 5(a) and (b)）, Such data would more clearly illustrate how the PLDOS varies spatially and spectrally near the BIC frequency.", "reviewer_3": "The authors present a near-field THz-TDs that enables the direct spatial mapping of BICs in metasurfaces. The results provide the evidence of a large partial LDOS (PLDOS) enhancement associated with BICs. In my view, this work presents several intriguing aspects that are valuable to readers and colleagues in this field. However, there are several points that require clarification.\n1. It should clearly explain the “nonlocal metasurfaces” to ensure that readers understand the nonlocal in this manuscript.\n2. The authors state that this is the first direct experimental evidence of LDOS enhancement by a BIC. However, in Ref. 24, an “optically inaccessible” BIC mode has been probed by electron energy-loss spectroscopy (EELS). They should compare this work to Ref. 24 in detail.\n3. The leakage of BIC states is generally observed by the quasi-BIC patterns. The author should more clearly present the influence of the probe on the BIC/quasi-BIC mode.\n4. In Fig. 4, the authors severally present calculated and measured spatial maps of PLDOS enhancement. Are these maps all normalized to the same maximum?\n5. The authors specifically explore the enhancement of PLDOS depending on different numbers of dimers N. whether the other parameters, for example, the thickness of gold film, will affect the enhancement?"}, "author_response": "We sincerely thank Reviewer #1 for the positive assessment of our manuscript and for recognizing the originality and significance of the experimental work. In both the manuscript and this response letter, we have clarified all the points raised by the reviewer. We believe these revisions substantially improve the manuscript’s clarity and impact.\n1. The manuscript states, \"we have developed a near-field THz time-domain spectrometer…\" suggesting the development of a new setup. However, unless I am mistaken, a rather similar system has been previously reported by the authors (e.g., ACS Photonics, 2022, 9, 5, 1547–1554). It would be more accurate to acknowledge that this established setup is now being utilized (or modified) to probe BICs in metasurfaces.\nResponse：We thank Reviewer #1 for pointing out and drawing the attention to our earlier publication (ACS Photonics 2023 10 (8), 2980-2986), in which we successfully demonstrated a near-field THz time-domain spectrometer for the quantitative determination of the local density of optical states (LDOS) in simple photonics systems, including planar interfaces made of gold, InSb, and quartz and also a resonant gold rod.\nBased on that established platform, we have demonstrated in this manuscript the spatially-resolved PLDOS of nonradiative and symmetry-protected BICs in the metasurfaces, which are invisible from the far-field measurement techniques. To clarify this distinction, we have revised the manuscript by adding:\n“Building on our previously developed dual near-field photoconductive-antenna-based THz time-domain microscopy 1, we have applied it to probe the PLDOS of BICs in metasurfaces that are not accessible from the far-field.” (Page 3 in the revised manuscript)\nAnd\n‘To directly measure BICs, we employ two THz microprobes to locally excite and probe the near-field response of the metasurface, where the microprobes function as sub-wavelength photoconductive antennas, as demonstrated in Ref. [20].’ (Page 5 in the revised manuscript)\n2. The concept of PLDOS is central to this study. It first appears in the abstract without explanation (whereas the LDOS is introduced well). Given its nuanced nature, a clear definition accompanied by relevant literature would greatly benefit readers. Unless I am mistaken, the current manuscript lacks any references on PLDOS, which should be addressed. For example, I note: Ly, Ousmane, et al. \"Partial local density of states from scanning gate microscopy.\" Physical Review B, vol. 73, no. 23, 2006, p. 235305. DOI: http://doi.org/10.1103/PhysRevB.73.235305.\nResponse: We thank Reviewer #1 for highlighting the importance of clearly defining PLDOS in the context of our study. We have to admit that our previous manuscript lacks a more explicit introduction to PLDOS, particularly given its central role in our work. We have added the following clear definition of the PLDOS in the introduction section of the revised manuscript, emphasizing that it represents the projection of the LDOS along the direction of the dipole moment of the source:\n‘Importantly, the dipole orientation strongly determines the spatial variations in its decay rate, especially in resonant subwavelength systems. Dipoles with different orientations interact with distinct combinations of polarized electromagnetic fields. Therefore, the partial (or projected) local density of optical states (PLDOS) becomes the more relevant quantity to describe light-matter interactions [1]. For instance, a dipole aligned along a unit vector 𝒆 interacts with 𝒅 modes whose electric fields are aligned in the same direction. At this level of detail, the PLDOS, denoted as 𝜌 (𝒆 ,𝒓 ,𝜔) , weighs each mode by the electric field projected along the dipole 𝑝 𝒅 𝟎 direction 𝒆 . The PLDOS captures the directional dependence of the LDOS, enabling the 𝒅 decomposition into contributions from specific polarizations and spatial orientations..’ (Page 2 in the revised manuscript).\nIn addition, we have now included a relevant manuscript as a reference:\n[1] W.L. Barnes, S.A.R. Horsley, W.L Vos, “Classical antennas, quantum emitters, and densities of optical states” Journal of Optics 22 (7), 073501 (2020).\nThis manuscript is fully focused on the description of the LDOS and PLDOS in optical systems.\nTherefore, we find this reference better suited than Ly, Ousmane, et al. \"Partial local density of states from scanning gate microscopy.\" Physical Review B, vol. 73, no. 23, 2006, p. 235305, which describes the partial local density of electronic states.\n3. The manuscript mentions “a fixed 10 μm distance between the emitter and detector, as shown in Figure 2(a)”. However, this distance is not clearly depicted in the figure. In my experience, determining and maintaining such spacing in near-field experiments can be challenging. On a related note, I see that experimental data shown in Figure 5c appears to originate at 0μm. Could the authors elaborate on how this distance was measured or verified, including any associated uncertainties?\nResponse: We thank Reviewer #1 for pointing out the missing information regarding the determination of the spacing between the two probes and their height above the sample surface during the measurement. To clarify this point, we have added a very detailed description of the probe alignment procedure to the revised supplementary material of the manuscript. We refer to this supplementary material by adding the following text to the revised manuscript:\n“The details on alignment of microprobes can be found in the Supplementary Information (Figures S2-3).”(Page 5 in the revised manuscript)\nThe added supplementary material is the following:\nFigure R1 shows photographs of the setup. As shown in Figure R1a, there are two 3D motorized linear stages for mounting the sample and the micro-probe for emission, respectively. In addition, the detector is mounted on a 1D linear stage that moves along the z-axis. To ensure both probes are at the same height, we first use the method we refer to as “Approach -Touch - Release.” In this approach, the sample is slowly moved toward the detector probe. Contact is detected when the probe visibly touches the sample surface, as observed with an optical camera through the reflected light of the optical pump on the probe (Figure R1b). At the moment of contact, the movement of the stage is immediately halted and then retracted in 1 μm steps until the probe returns to its original shape. This position is defined as the zero-height reference point (i.e., 0 μm above the sample surface). The same process is repeated for the emitter to confirm its height over the metasurface.\nFigure R1 (a) Photo of the dual-photoconductive antenna-based THz near-field microscopy. (b) Camera view of the sample loading process: (i) Approach: moving the metasurface toward the microprobes with a motorized 3D linear stages; (ii) Touch: Sample touching the probe, which is observed by the reflection of the laser pump for THz generation/detection on the probe; (iii) Release: Movement backward of the probe by 10 µm.\nThe distance between the two probes is precisely controlled using the 3D motorized stage for the emitter. The THz signal from the detector is recorded as the emitter is scanned along the y-axis, from a position above to another position below the detector. This scanning allows us to identify the position of maximum signal intensity, which we define as the lateral center of the detector. To more accurately determine this central y-position, Figure R2 presents a spatial mapping of the emitter as it also moves over a 100 μm range along the x-axis. The results clearly show that for all x-positions, the maximum signal consistently occurs at approximately 138 μm along the y-axis (with the initial position labeled as ‘0’). This confirms the central position of the detector along the y-direction. Lastly, the emitter is moved along the x-axis toward the detector until physical contact is observed via the camera. At this point, the laser beam on the microprobes begins to shake, as shown in Figure R1b (ii), indicating the contact. This point is defined as the zero position along the x-axis. During the measurements, the emitter is displaced by 10 μm from the detector along the x- and y-axes to prevent the excitation laser of the emitter and reflector from undergoing multiple reflections between the closely spaced microprobes.\nFigure R2. (a) Illustration of the emitter movement over the xy plane to confirm the relative distance between x and y. (b) Spatial mapping of THz amplitude at 0.3 THz measured by the detector when the emitter is moved over xy plane. Here ‘0’ indicates the initial starting position for the emitter.\n4. Simulation results are placed in the supplementary material and seem to exhibit narrower features compared to experimental data. This discrepancy should be elaborated upon. Where possible, I would recommend to:\n(1) Integrate key simulation results into the main figures – particularly Figure 1b, where the most important feature takes up only a small portion of the figure.\nResponse: Thank you for the helpful suggestion. Following the recommendation of Reviewer #1, we have integrated the key simulation results into the main figures of the revised manuscript. Specifically, we have updated Figure 1 to include the simulated transmission spectra and Q-factor as a function of the incidence angle. In the revised manuscript, we have added the following description:\n‘The off-normal excitation breaks the inversion symmetry of the system, causing the BIC to evolve into a quasi-BIC that becomes accessible from the far field. As the angle of incidence increases, the resonance redshifts and broadens, indicating a reduction in Q-factor. These spectral changes are more clearly observed in the relative transmission spectra, normalized by the transmission at normal incidence, as shown in the right panel of Figure 1(b). The observed redshift and Q-factor reduction are qualitatively consistent with simulation results in Figure 1(c), obtained using the scuff-transmission method in SCUFF-EM (see Methods). In these simulations, the BIC in an infinite array of rod dimers appears at 0.395 THz and exhibits an infinite Q-factor. We note a slight mismatch between the experimental and simulated resonance frequencies. This mismatch arises from idealized assumptions in the simulations with a lossless, infinite, and homogeneous environment, that differ from real experimental conditions (a more detailed discussion is provided in the Supplementary Information). Despite this discrepancy, the BIC framework remains an effective tool for the design of nonlocal metasurfaces.’ (Page 5 in the revised manuscript)\nFigure R3. Sample geometry and far-field transmission. (a) Illustration of a unit cell of the metasurface and transmission measurements with changing incidence angles. The inset shows an optical microscope image of a unit cell of the fabricated metasurface, where the separation between the rods is d = 80 μm. (b) The left panel displays the THz transmission amplitude of the metasurface for incidence angles of 0°, 20°, 30°, and 40°. The right panel shows the corresponding relative transmission (Rt) spectra to the transmission at normal incidence. These normalized spectra highlight the emergence of a quasi-BIC mode, which becomes prominent under off-normal incidence. In these measurements, the THz electric field is polarized along the y-direction (s-polarized). (c) Calculated transmission spectra as a function of the incidence angle. The symbols and solid line represent the extracted Q-factor of the quasi-BIC mode by 𝑄 =𝑓⁄∆𝑓, where 𝑓 is the resonant frequency of the quasi- 0 0 BIC mode, and ∆𝑓 represents the linewidth (half-width at half-maximum) of the quasi-BIC mode in the Rt spectra retrieved via Fano fitting.\n(2) Overlay experimental and simulation data on the same axes for direct comparison (e.g., in Figure 3d), as independently scaled axes and the separation into insets may obscure important differences.\nResponse: Thank you for pointing this out. In the revised manuscript, we have overlaid the experimental and simulation data used the same axis scale for different graphs to enable direct comparison, as shown in Figure R4. This revision provides a clearer visualization of the similarities and differences between the two datasets.\nAs illustrated in the updated figure, both experimental and simulated results consistently show that the BIC mode becomes sharper with an increasing number of unit cells. This behavior arises from the collective response of a larger number of rods. However, we acknowledge that, due to the limited frequency resolution of our measurement system, the experimentally observed peaks appear broader compared to those obtained from simulations.\nIn the revised manuscript, we have updated Figure 3 as the one shown below,\nFigure R4. Enhancement of the PLDOS by the BIC. (a-c) Calculated and measured PLDOS enhancement 𝑃 /𝑃 as a function of frequency and position along the y-direction at half of the width of one rod of the unit 𝑦𝑦 0 cell, indicated by the dotted line in the inset of (a). The calculated PLDOS enhancement for a single dimer (a), and for an array consisting of 9 × 9 dimers (b). (c) Measured PLDOS enhancement of the metasurface consisting of 84 × 84 unit cells. (d) Measured and calculated frequency-dependent PLDOS enhancement at the edge of a rod (dashed-horizontal lines marked in (a-c)) in arrays with different numbers of dimers 𝑁×𝑁: (i) 84 × 84 (measurement only), (ii) 9 × 9, (iii) 7 × 7, (iv) 5 × 5, (v) 3 × 3, and (vi) 1 × 1. (ii-vi) are calculations. The insets in (iv-vi) are measurements of the finite-size arrays. The peaks in the PLDOS enhancement for dimer arrays (N > 1) result from the collective nature of the BIC confined in the finite array. (e) Spatial distribution of the PLDOS enhancement for the 9×9 array at the respective PLDOS peaks, corresponding to the labeled peaks in the inset figure d, panel (ii). (f) Calculated quality factor Q, and (g) PLDOS enhancement extracted from the lowest frequency peak as a function of the number of dimers in the array (averaged over 10 positions near the rod edge, i.e., over 10 µm along the y-direction). The solid lines represent exponential fits to the data.\n(3) Apply color scales that do not obscure data, for example, Figure 4 appears to contain saturated / truncated colourmap data, and seems to have been interpolated instead of showing the actual data in each measurement pixel.\nResponse: Thank you for the valuable comment. In the revised manuscript, we have adjusted the color scales in Figure 4 to avoid saturation and truncation, ensuring that the full dynamic range of the data is visible. Additionally, we have updated Figure 4 to display the actual measured data at each pixel without interpolation, providing a more accurate representation of the experimental results.\n(4) Each subfigure in Figure S1 b, c, have different frequency axes.\nResponse: In the revised manuscript, Figure S1c from the supplementary materials has been moved to the main text as Figure 1c to highlight the incidence angle-dependent quasi-BIC mode more clearly. To facilitate direct comparison between the simulation and experimental results, we have used the same frequency scale in Figure 1c and the inset of Figure 1b. Compared to the experimental data, the simulated quasi-BIC mode exhibits a sharper resonance dip and a larger redshift with increasing incidence angle. This discrepancy arises because, in the simulations, the gold rods were modeled as perfectly planar rectangles embedded in a uniform dielectric medium with a refractive index of n = 1.55, and material losses were neglected. The value n = 1.55 was chosen as an average of the refractive indices of air (n = 1.0) and the quartz substrate (n = 2.1).\nTo validate this explanation, we have added new simulation results in the revised supplementary information, as shown in Figure R5 (Figure S1 in the Supplementary). Figure R5a presents the measured transmission spectra under plane-wave excitation with varying incidence angles from 0° to 40°. To enhance the visibility of the quasi-BIC mode, all oblique incidence spectra are normalized by the normal incidence spectrum, and presented as the relative transmission (Rt) spectra. Figure R5b compares the Rt spectra from the experiment with two sets of simulation results: one assuming the perfect electric conductor embedded in a homogeneous medium, and the other modeling the rods at the interface between air and quartz, with material losses included. The latter scenario more accurately reflects the experimental conditions. It is clearly shown that the frequency of the quasi-BIC mode is redshifted with increasing incidence angle and its linewidth gets broader. These simulations show that the frequencies of the quasi-BIC modes align much more closely with the experimental results, confirming the significant influence of both the dielectric interface and the material losses. Despite differences in the dispersion behavior between the two simulation approaches, all cases consistently indicate the emergence of the BIC mode near 0.395 THz.\nFigure R5. Incidence angle-dependent transmission spectra of the metasurface. (a) Normalized transmission spectra of the metasurface under different illumination angles, from 0° to 40° for y-polarized radiation. (b) Relative transmission spectra of the metasurface under oblique incidence normalized by the transmission under normal incidence: (i) experiment, (ii) simulation of infinite size metasurface of perfect electric conductor (PEC) in an homogeneous environment, (iii) simulation of infinite metasurface of gold rods on a quartz substrate.\n(5) A version of Fig. S1c using experimental and the same frequency axes would be beneficial (the data is already there.) Alternatively, a version of Fig. S1b with simulation data might help comparisons (but is less effective at capturing the emergence of coupling to the BIC mode).\nResponse: Thank you for the good suggestions. We have modified Fig. S1 in the revised manuscript, which is shown in Figure R5 above. Following your recommendation, we have revised Figure S1 in the revised manuscript to enable a direct comparison, as shown above. Specifically, in the new version of Fig. S1b:\nPanel (i) shows the experimental relative transmission (Rt) spectra.\nPanel (ii) presents the simulation results under idealized conditions, where gold rods are modeled as perfectly rectangular structures embedded in a homogeneous dielectric medium (n = 1.55), with material losses neglected.\nPanel (iii) shows simulations that more closely reflect experimental conditions by including material losses and modeling the rods at the air–quartz interface.\n5. The study employs a mapping step size of 4 μm, while to my knowledge the photoconductive antennas typically have an effective resolution of 10 μm according to the specifications. Could this lead to oversampling and feature blurring? Do the authors have an estimate of the system's point spread function (PSF), or can they comment on the effective resolution?\nResponse: We thank Reviewer #1 for the insightful comment regarding the mapping step size and effective resolution of the system. While it is true that photoconductive antennas (PCAs) typically have an effective resolution on the order of ~10 μm due to the carrier diffusion in the substrate, we chose a step size of 4 μm to ensure accurate sampling of the signal. Although this results in oversampling relative to the nominal resolution, it does not lead to feature blurring; rather, it improves the robustness of spatial mapping and enables more accurate numerical post- Figure R6. (a) Optical microscope image of metallic strips with different widths from 5 μm to 35 μm with a step of 5 μm. During the measurement, the microprobes are placed at a height of 10 μm over the sample surface. (b) THz near-field intensity spectra in arbitrary units measured in a line scan over the metal strips. (c) THz field intensity of the line scanning over the metal strips at 0.88 THz. The yellow regions represent the metal strips. processing, such as interpolation or gradient-based analysis.\nTo explicitly address the effective resolution of the setup in the revised manuscript, we have determined the line spread function (LSF). The LSF characterizes how an optical system responds to a narrow, line-shaped input and serves as the one-dimensional analogue of the PSF. We have fabricated a line array of metal strips with widths ranging from 5 μm to 35 μm and a thickness of 100 nm, as shown in Figure R6(a). Each strip is 1 mm in length. A line scan was performed across the array using the THz near-field microscope with a step size of 5 μm and the probe at a height of 10 μm. The lateral separation between the two microprobes was 20 μm along the direction of the line scan (x) and 10 μm perpendicular to this direction (y). Figure R6(b) presents the measured THz intensity spectra at each scan position. This field enhancement remains clearly detectable even for the 5 μm-wide metal strips, highlighting the system’s ability to detect subwavelength metal features, as shown in Figures R6(b) and R6(c).\nThe linewidth of the intensity peak obtained by fitting the averaged THz intensity profile (integrated between 0.7 THz and 1.0 THz) with a Gaussian function in Figure R6(c), appears broader than the intrinsic spatial resolution of the microprobes. This broadening is primarily attributed to the ~20 μm lateral separation between the two microprobes used during the measurement.\nThe spatial resolution description of the system has been added to the supplementary information in the revised version of the manuscript. In the maintext, we have added the following description of the spatial resolution for THz near-field microscopy:\n‘While the line spread function of the near-field probes has been measured to be 20 μm (see Supplementary Information, Figure S4), we choose a smaller step size for a detailed mapping’ (Page 6 in the revised manuscript)\n6. Figure 3 illustrates the trend of increasing BIC behavior with metasurface size. However, experimental data for N = 7 and N = 9 are absent. Have measurements been conducted for these configurations? Additionally, could signal processing techniques, such as zero-padding, enhance the visibility of spectral features, perhaps in smaller metasurfaces where resonances are wider (e.g., N = 5)?\nResponse: Thank you for the questions. Measurements for N = 7 and N = 9 have indeed been conducted and are now included in Figure 3 of the revised manuscript. In the previous version of the manuscript, these data were quite noisy, with multiple peaks due to echoes from the thin substrate of 500 μm thickness. To solve this issue, we have prepared new samples on a thicker substrate, with a thickness of 5 mm, and repeated all the measurements on the size-dependent metasurfaces.\nFigure R7. Measured PLDOS enhancement of the metasurfaces with finite size from N = 1 to N = 9.\nFor a more detailed presentation of the measured PLDOS spectra across different metasurface sizes, please refer to Figure R7. The measured value of the PLDOS enhancement is not quantitatively comparable with the calculated value, one reason is that the simulations assume ideal, lossless materials, whereas material losses in the actual samples contribute to peak broadening and the reduction of the PLDOS in the measurements. Another reason is due to the frequency resolution in our setup, which is approximately 20 GHz, determined by a 50 ps time window of the measurements. The sharp resonance peaks predicted by the simulations for N = 7 and N = 9 can not be resolved in the experimental data. Despite these limitations, a clear trend is evident in the experimental results for metasurfaces with increasing size. As the number of unit cells increases, the BIC mode becomes progressively sharper, and the PLDOS enhancement becomes more pronounced—consistent with the expected collective behavior predicted by theory and simulation.\nFurthermore, we thank Reviewer #1 for the suggestion to use zero-padding to improve the signal. We have applied zero-padding as a signal processing technique to enhance the spectral resolution. We compare different levels of zero-padding, varying the number S of zeros appended to the data, in Figure R8. Although zero-padding helps smoothing the spectra and slightly improves the apparent spectral resolution, it is inadequate for resolving the narrow resonances observed in the simulations. This limitation stems from the finite time window of our raw data, which consists of only 1000 points sampled at 0.05 ps intervals, corresponding to a total duration of 50 ps, and thus setting a limit on the achievable frequency resolution.\nIn the revised manuscript, we have added a description about the applied zero-padding to the transients to enhance the visibility of spectral features:\n‘Owing to the limited scan range of the time-delay stage, zero padding is employed to the time-domain data prior to Fourier transformation, for enhancing the frequency spectral resolution in the frequency domain.’ (Page 15 in the revised manuscript)\nFigure R8. Comparison of PLDOS enhancement spectra data without zero-padding and with varying the number of zeros appended to the THz transients.\n7. Are there simulations corresponding to Q and Pyy/Po in Fig. 3(f,g)?\nResponse: Thank you for the question. In fact, the results presented in Figure 3(f, g) are simulations to better illustrate the BIC’s properties with infinite Q-factor and the saturation of the PLDOS enhancement with the metasurface’s size. These simulations result in values of Q and Pyy/Po that are larger than those retrieved from the experiments as shown in Figure R9, due to the limitation of the frequency resolution (20 GHz) in the spectra, the materials’ loss and the effect of substrate, which are neglected in simulations.\nFigure R9. The extracted Q-factors and PLDOS enhancement of the BIC mode from the measurement of the metasurfaces with finite sizes composed of different numbers of dimers.\nOn page 15 of the revised manuscript, we have added the following detailed description of the experimental setup, including the scanning range of the time-delay stage, which determines the time window of the THz transients and the resulting frequency resolution of the spectra in the THz near-field microscopy system. The measured results have been added to the Supplementary Information (Figure S9).\n8. I suggest making a small note clarifying that quartz is used as a reference due to it being the substrate material.\nResponse: We thank Reviewer #1 for this helpful suggestion. As noted in the manuscript (page 7), we have already stated that “The reference sample that is used for normalizing the measurements of the array is a quartz substrate.” This choice is indeed motivated by quartz being the substrate material. We have clarified this point explicitly in the revised manuscript to avoid any ambiguity by adding the following sentence to the revised manuscript:\n‘A bare quartz substrate, identical to the sample substrate, was used as the reference to normalize the measurements and isolate the optical response of the metasurface’. (Page 5 in the revised manuscript)\n9. Further comments on potential related literature:\nThis work sits within a broader research field using near-field antennas to probe photonic structures (wherein the present paper considers the case where two antennas are very close together). The authors are therefore encouraged to consider whether any of the following works could be relevant and might merit citation:\n- Gansch et al., Light: Sci. Appl., 2016. DOI: https://doi.org/10.1038/lsa.2016.147\n- Yermakov et al., Phys. Rev. B, 2018. DOI: https://doi.org/10.1103/PhysRevB.98.195404\n- Bulu et al., New J. Phys., 2005. DOI: https://doi.org/10.1088/1367-2630/7/1/223\n- Li et al., Phys. Rev. Applied, 2025. DOI: https://doi.org/10.1103/PhysRevApplied.23.014084\n- Lezhennikova et al., Comm. Physics 2023, https://doi.org/10.1038/s42005-023-01347-1\n- Lezhennikova et al., ACS Photonics 2023, https://doi.org/10.1021/acsphotonics.3c01312\nAt this point I would like to declare that I am an author of the last paper of the above list (Lezhennikova et al., ACS Photonics 2023). It presents a related approach using the same combination of antennas (albeit at different locations on the opposite side of a photonic crystal sample), and it may offer complementary insights, particularly regarding a spatial Fourier-space analysis and the observation of enhancement / suppression of a near terahertz field emitter as a result of the available photonic states. Of course, I am not requesting that this work be citated in the revised manuscript, but merely making the authors aware of it in case they happen to consider it relevant or useful.\nResponse: We sincerely thank the reviewer for highlighting these relevant contributions and for providing valuable context regarding our study's position within the broader field of near-field antenna probing of photonic structures. We found that several of the suggested papers offer useful context and complementary insights into our approach. In particular, we have cited the works by Gansch et al.( Light: Sci. Appl., 2016. DOI: https://doi.org/10.1038/lsa.2016.147), Lezhennikova et al. (Comm. Physics 2023, https://doi.org/10.1038/s42005-023-01347-1, ACS Photonics 2023, https://doi.org/10.1021/acsphotonics.3c01312), and Li et al (Phys. Rev. Applied, 2025. https://doi.org/10.1103/PhysRevApplied.23.014084) in the revised manuscript.\nWe also appreciate that the Reviewer has made us aware of the study by Lezhennikova et al. (ACS Photonics, 2023), which provides an example of near-field terahertz spectroscopy applied to photonic crystal slabs. This work demonstrates emission enhancement and suppression linked to the photonic band structure, and it is especially relevant. Both our study and that of Lezhennikova et al. employ broadband terahertz dipole antennas in a near-field geometry to investigate structured photonic systems. While our setup uses co-located emitter and detector antennas with subwavelength separation on the same side of the sample to directly probe localized field enhancements and PLDOS enhancement, Lezhennikova et al. position the antennas on opposite sides to map transmitted near-field patterns and directional emission behavior. These distinct configurations access complementary physical observables. Despite the methodological differences, both approaches highlight the utility of near-field techniques in accessing photonic phenomena not readily observable in the far field. Together, they underscore the richness of structured photonic environments and the diverse ways in which near-field probes can be used to explore them.\nIn the revised manuscript, we have added the following description in the introduction,\n‘With the advancement of near-field techniques spanning from the visible to the terahertz (THz) regime, these methods have become widely utilized for controlling and probing electromagnetic waves in engineered photonic and metamaterial platforms [27–30]. Among them, THz near-field microscopy integrated with dual microprobes offers a powerful approach for investigating diverse photonic phenomena, including spatial coherence and Fourier components of resonant modes [25, 31]’ (Page 3 in the revised manuscript )\nWe have included the following references in the introduction of the revised manuscript:\n(1) Gansch et al., Light: Sci. Appl., 2016. https://doi.org/10.1038/lsa.2016.147\n(2) Lezhennikova et al., Comm. Physics 2023, https://doi.org/10.1038/s42005-023-01347-1\n(3) Lezhennikova et al., ACS Photonics 2023, https://doi.org/10.1021/acsphotonics.3c01312\n(4) Li et al., Phys. Rev. Applied, 2025. https://doi.org/10.1103/PhysRevApplied.23.014084\n(5) Bohn B et al., Nano Letters, 2015 https://doi.org/10.1021/acs.nanolett.5b00692\n10. I think it would be broaden the appeal of this work to speculate more, in the conclusions, about the implications of this work at other frequency ranges, or for photonic- and quantum- technologies in general, e.g., in the context of quantum or single photon emitters (if appropriate), and beyond, at the authors’ discretion.\nResponse: Thank you for the suggestion. We agree that broadening the scope of the discussion can help highlight the broader impact of our work. In the revised manuscript, we have expanded the conclusions to include speculative remarks on the potential implications of our findings across other frequency ranges and emerging technologies.\nThe following text has been added in the conclusions of the revised manuscript:\n’Owing to the scale invariance of Maxwell’s equations, the design principles and insights developed here for terahertz metasurfaces can be naturally extended to the infrared and visible regimes. In particular, the ability to engineer BIC modes and enhance the PLDOS opens up promising avenues for controlling light–matter interactions in nanophotonic systems. This has important implications for photonic and quantum technologies, such as the enhancement of the emission efficiency of single-photon sources, the tailoring of the spontaneous emission rates in cavity quantum electrodynamics, or the improvement of nonlinear optical responses in integrated photonic platforms.’ (Page 13 in the revised manuscript )\n11. Minor Suggestions:\n- The measurement conditions (“4 μm resolution, 10 μm height”) are repeated more than once.\nResponse: The repeated explanation of the experimental method and simulation settings have been removed in the revised manuscript.\n- Clarify the method for Q-factor extraction related to the phrase “The Q-factor was extracted from the first peak (fundamental resonance) of the frequency-dependent PLDOS enhancement”\nResponse: Thank you for your helpful reminder, and we apologize for the omission in the original explanation regarding the Q-factor extraction.\nWe extracted the Q-factor from the frequency-dependent PLDOS enhancement spectra, taking the metasurface with finite size 𝑁 = 5 as an example (see Figure R10). The first peak in the spectrum is identified as the fundamental BIC mode, while the higher-order peaks are attributed to collective hybridized modes formed through the coherent coupling of localized resonances across the unit cell array. These resonances increase in number and collapse into the BIC mode when the size of the array increases.\nTo extract the Q-factor, we performed a Fano line-shape fitting on the first resonant peak using the following equation in Eq.1,\n2\n(𝑞+\n2(𝑓−𝑓0)\n)\n𝑦 = 𝐴× ∆𝑓 +𝐵, (1)\n2\n1+(\n2(𝑓−𝑓0)\n)\n∆𝑓\nwhere 𝑓 is the resonant frequency of the peak, ∆𝑓 is full width at half maximum (FWHM), 0 and 𝑞 is Fano asymmetry parameter. In addition, A and B are scaling and offset constants, respectively. The quality factor is then calculated as: 𝑄 = 𝑓 /∆𝑓 . This method provides an 0 accurate characterization of the fundamental resonance associated with the BIC mode.\nThis explanation has been included in the supplementary of the revised manuscript.\nFigure R10 Extracted Q-factor obtained by fitting the first resonant peak in the PLDOS enhancement spectrum using a Fano line profile. The example shown corresponds to a metasurface of finite size 𝑁=5.\n-“We have calculated and the Q-factor…” – is this a typo?\nResponse: We are sorry for the typo. It has been corrected in the revised manuscript. We have further checked the manuscript, and other typos have been corrected.\n- Please specify the model of the Protemics GmbH antenna used, as there are several and they vary in characteristics.\nResponse: The antenna model we used for this measurement is the TeraSpike TD-800-X-HR-WT. This information has been added in the methods of the revised manuscript.\nIn conclusion, this study presents impressive results in a technically challenging field. The observations are significant, and the technique holds promise. My suggestions aim to improve clarity, transparency, and context. Congratulations to the authors on the experimental achievement.\nResponse: We sincerely thank the reviewer for the thoughtful and encouraging comments. We greatly appreciate the recognition of the technical challenges involved in our work, as well as the significance of our observations and the potential of our technique. We also value the reviewer’s constructive suggestions, which we have carefully addressed in the revised manuscript to improve its clarity and transparency.\nReviewer #2 (Remarks to the Author):\nIn the manuscript “Near-field Probing of the Local Density of Optical States Enhanced by Bound States in the Continuum in Nonlocal Metasurfaces” by Jie Ji et al, the authors demonstrate the near-field mapping of BICs in metasurface and vertify the associated PLDOS enhancement. Unfortunately, I find that this manuscript demonstrates some overlaps with the prior works (ter Huurne, Stan, et al. \"Direct Observation of Lateral Field Confinement in Symmetry‐Protected THz Bound States in the Continuum.\" Advanced Optical Material 11.6 (2023): 2202403; van Hoof, Niels JJ, et al. \"Unveiling the symmetry protection of bound states in the continuum with terahertz near-field imaging.\" ACS photonics 8.10 (2021): 3010-3016.). As a result, the results in this manuscript are not particularly innovative and impressive. Regrettably, I cannot recommend this manuscript for publication in Nature Communications.\nResponse: We sincerely thank the reviewer for the assessment. However, we respectfully disagree with the statement that prior work in a particular system removes the novelty of new work. Any field in science is plagued with high-impact results in the same or similar systems. Future breakthroughs are built upon previous results. The distinct novelty of this work, as compared to previous references, can be summarized by:\n First direct PLDOS mapping of BICs: While the prior works (van Hoof et al., ACS Photonics 2021; ter Huurne et al., Adv. Optical Mater. 2023) provided pioneering near-field imaging of symmetry-protected BICs at THz frequencies and characterized lateral confinement and lifetime of BIC modes, our manuscript presents, to the best of our knowledge, the first direct experimental mapping of the partial local density of optical states (PLDOS) enhanced by BICs. This is a critical advancement since LDOS (and PLDOS) directly quantifies the potential for enhanced light–matter interaction, a parameter not measured or addressed in the prior studies.\n Observation of the saturation of the PLDOS despite the diverging Q: Our findings show that, although the quality factor of the BIC diverges with increasing metasurface size, the PLDOS reaches a finite limit. This nuanced result, not discussed in the literature, has significant implications for the theoretical limits of BIC-enhanced devices, demonstrating that mode volume and nonlocal interaction fundamentally cap the achievable enhancement.\n Distinctive Focus on PLDOS Physics and Application Relevance: The prior THz works focused on lateral field confinement, symmetry-breaking behavior, and near-field excitation schemes in metallic resonator arrays. In contrast, our work establishes a direct link between symmetry-protected BICs and enhanced LDOS, which is a critical metric for photonic engineering, quantum optics, and sensing applications. This shift in focus represents a substantive conceptual innovation.\nTo further stress the novelty of our results, we have added the following text to the introduction of the revised manuscript in page 3:\n‘Due to these limitations, LDOS enhancements are commonly indirectly obtained from measurements of the lifetime or decay rate of excited states [21–23]. However, this approach is not applicable to BICs, as they do not couple to the far-field continuum. Instead, near-field techniques with specific excitation or detection methods, are required to investigate these dark modes and visualize their spatial field distribution [24–26]. or example, Dong et al. used electron energy loss spectroscopy to probe the near-field response of optically inaccessible BIC modes in silicon nano-antenna arrays, revealing strong spatial confinement of the associated electromagnetic energy at nanometer resolution [24]. Hoof et al. used a double THz near-field probe microscope to measure the near field amplitude and phase of BICs and reveal the symmetry protection that suppresses radiation losses [25] With the advancement of near-field techniques spanning from the visible to the terahertz (THz) regime, these methods have become widely utilized for controlling and probing electro- magnetic waves in engineered photonic and metamaterial platforms [27–30]. Among them, THz near-field microscopy integrated with dual microprobes offers a powerful approach for investigating diverse photonic phenomena, including spatial coherence and Fourier components of resonant modes [25, 31, 32].\nIn this work, we present the first direct experimental evidence of the enhanced PLDOS associated with BICs in metasurfaces. Building on our previously developed dual near-field photoconductive-antenna-based THz time-domain microscopy [20], we have applied it to probe the PLDOS of BICs in metasurfaces that are not accessible from the far-field.’ (Page 3 in the revised manuscript )\nWe hope this clarifies how our current work builds upon—but goes significantly beyond—previous efforts in this area. Thank you again for your consideration.\nAdditional comments are as follows:\n1. In the work, the authors mentioned, “To investigate the far-field properties of the metasurface and validate the presence of the BIC, we measured the amplitude transmission spectra T(ω, θ) through the sample at various angles of incidence θ with a far-field THz-TDS system.” They demonstrate the existence of BIC through angle-resolved far-field experiments. However, the limited comparison of experimental results under only two angular incidence is not sufficient. More importantly, rigorous validation criteria for BIC require explicit verification of the presence of polarization topological charges.\nResponse: We thank the reviewer for their valuable observations. Rigorous validation of the investigated BIC was provided in the manuscript “Spectral and temporal evidence of robust photonic bound states in the continuum on terahertz metasurfaces” Optica 6, 996 (2019). Using an electric/magnetic coupled dipole model, it was demonstrated in this manuscript the suppression of the imaginary component of the eigenvalue of the metasurface at the BIC frequency when the two rods are equal in length, which is a ubiquitous indication of an infinite lifetime in a lossless system, thus a BIC. To clarify this point, we have added the following text on the revised manuscript:\n‘Given the designed dimensions, we expect a BIC mode emerging at a frequency of 0.395 THz as demonstrated in Ref. [35]. This is supported by a theoretical analysis using a coupled detuned-dipole model, which shows that the imaginary component of the system’s eigenvalue vanishes at this specific frequency.’ (Page 4 in the revised manuscript )\nRegarding the angle-resolved far-field measurements: although only two incidence angles are presented in the main text to illustrate the emergence of the quasi-BIC, we have performed a comprehensive set of measurements across a broad range of incidence angles. These results, included in the Supplementary Information (Figure S1) in the previous version of the manuscript, clearly demonstrate the evolution of the quasi-BIC resonance with angle, showing its progressive narrowing toward normal incidence and its disappearance at θ = 0°, consistent with the behavior of symmetry-protected BICs. To make this more evident in the main text, we have added two additional angle-dependent measurements in Figure R1, as well as, incidence-angle-dependent simulation results as shown in Figure R2.\nFigure R1. Sample configuration and far-field transmission. (a) Illustration of a unit cell of the metasurface and transmission measurements with changing incidence angles. The inset shows an optical microscope image of a unit cell in the fabricated metasurface, where the separation between the rods is d = 80 μm. (b) The left panel displays the THz transmission amplitude of the metasurface for incidence angles of 0°, 20°, 30°, and 40°. The right panel shows the corresponding relative transmission (Rt) spectra to the transmission at normal incidence. These normalized spectra highlight the emergence of a quasi-BIC mode, which becomes prominent under off-normal incidence. In these measurements, the THz electric field is polarized along the y-direction (s-polarized). (c) Calculated transmission spectra as a function of the incidence angle. The symbols and solid line represent the extracted Q-factor of the quasi-BIC mode by 𝑄 =𝑓⁄∆𝑓, where 𝑓 is the resonant frequency of the quasi-BIC mode, and ∆𝑓 represents the linewidth (half-width at half-maximum) of the quasi-BIC mode in the Rt spectra via Fano fitting.\nFigure R2. Incidence angle-dependent transmission spectra of the metasurface. (a) Normalized transmission spectra of the metasurface under different illumination angles, from 0° to 40° for y-polarized radiation. (b) Relative transmission spectra of the metasurface under oblique incidence normalized by the transmission under normal incidence: (i) experiment, (ii) simulation of infinite size metasurface of perfect electric conductor (PEC) in an homogeneous environment, (iii) simulation of infinite metasurface of gold rods on a quartz substrate.\n2. The authors assert direct observation of a BIC mode through near-field measurements.\nWhile the microscopic probe as an excitation source in Figure 2(a) indicates a tilted orientation, which likely excites a quasi-BIC rather than a genuine symmetry-protected BIC. If the tilt angle of microprobe emitter is systematically varied, would the measured radiative leakage (e.g., resonance linewidth and Q-factor) show a significant variation—a hallmark of quasi-BIC behavior?\nResponse: We appreciate the reviewer’s thoughtful question and the opportunity to clarify the nature of near-field excitation in our experiment. In our study, the BIC mode is excited by a near-field point dipole source (photoconductive antenna) placed in close proximity (~10 μm) to the metasurface, as shown in Figure R3(a). Unlike plane-wave excitation, which is subject to well-defined incidence angles and symmetry constraints, a point dipole emitter inherently lacks translational symmetry and exhibits a broad angular spectrum. As a result, this localized source can couple to all the available modes, including those that are symmetry-protected from far-field radiation. Importantly, this localized near-field excitation does not convert the BIC into a quasi-BIC. Rather, it leverages local symmetry breaking at the excitation site, without altering the intrinsic radiation symmetry of the mode. Since the far-field radiation channel remains forbidden due to the global symmetry of the structure, the excited BIC mode continues to exhibit a high Q-factor, consistent with its bound nature.\nRegarding the orientation of the emitter, although the near-field dipole is slightly tilted (as shown in Figure R3(b), this tilt does not significantly affect the 𝐸𝑦 component of the emitted THz electric field, the component responsible for coupling to the BIC in the metasurface. Therefore, the tilt angle θ does not appreciably alter the coupling efficiency, resonance linewidth, or Q-factor of the BIC mode under this specific excitation geometry.\nTo clarify these points, we have added Figure R3 to the Supplementary Information of the revised manuscript, in addition to a section, explaining the local excitation of the BIC and the negligible effect of tilt angle.\nFigure R3. (a) Illustration of the THz near-field measurement on the metasurface. (b) The sketch of the direction of THz electric field direction of the emitter and detector corresponding to the rods during the measurement.\n3. In Figure 3(d), significant discrepancies are observed between the experimental and simulated curves, including variations in linewidths and resonance frequency shifts. Particularly for N = 5, anomalous multi-peaks are in the measured results. The authors should explain these differences in detail. Although the authors attribute peak broadening to the system’s frequency resolution limitations, they have not provided critical parameters about the delay window and the spectral resolution of their system. In addition, the experimental spectra in Figure 3(d) exhibit insufficient sampling points (the inflections in the experimental curve are overly pronounced). Since the signal fluctuation is also influenced by noise during measuring, this effect becomes particularly pronounced when the number of sampling points is limited.\nResponse: We thank Reviewer #2 for the thoughtful and detailed comments. We acknowledge the discrepancies between the experimental and simulated curves in Figure 3(d), particularly the broader linewidths, resonance frequency shifts, and the appearance of multiple peaks in the experimental data for N = 5, (𝑁×𝑁 being the total number of unit cells) These differences primarily arise from limitations in our experimental setup and practical constraints not captured in the idealized simulations. The contributing factors to the discrepancies are the following:\n(1) Spectral resolution: The frequency resolution of our terahertz time-domain spectroscopy (THz-TDS) system is approximately 20 GHz, determined by a time-domain delay window of 50 ps with a step size of 0.05 ps. This resolution is not sufficient to fully resolve the sharp spectral features predicted by simulations, particularly for metasurfaces with larger unit cell numbers (e.g., N = 7, and 9), where high-Q BIC modes are expected. Nevertheless, the trend of increasing PLDOS enhancement with metasurface size remains evident and consistent with simulations, as shown in Figure R4. Moreover, while simulations assume ideal, lossless materials, the actual samples experience intrinsic losses, which further broaden the resonances in the experimental results.\nFigure R4. The measured PLDOS enhancement of the metasurface with finite size from N = 1 to N = 9.\n(2) Anomalous Multi-Peaks for N = 5: The multiple peaks observed in the experimental data for N = 5 are likely caused by interference effects from echo signals reflected within the sample’s substrate, which has a thickness of 500 µm. To avoid this problem, we have prepared new samples using thicker substrates (5 mm thickness) and repeated the measurements of PLDOS enhancement in the finite-size structures. In the revised manuscript, all relevant data have been updated accordingly, as presented in Figure R4.\nWe have revised the manuscript to include these critical parameters in the methods section by adding the sentence:\n‘The frequency resolution of the terahertz near-field time-domain spectroscopy (THz-TDS) system is approximately 20 GHz, determined by the 50 ps time window from the delay stage.’ (Page 14 in the revised manuscript )\nIn addition, we have clarified the potential sources of discrepancy between experiment and simulation. Despite the experimental limitations, the overall trend remains clear: as the number of unit cells increases, the BIC mode becomes sharper, and the PLDOS enhancement more pronounced, consistent with the expected collective behavior.\n4. Regarding the Q factor variations in Figure 3(f), the authors claim that the Q factor is defined as the reciprocal of the mode’s lifetime. They should provide the detailed calculated method of the Q factor and display the experimental trend of Q-factor evolution with the number of dimers (N).\nResponse: We thank the reviewer for pointing out the need for clarification regarding the calculation of the Q factor and its dependence on the number of dimers 𝑁×𝑁. We extracted the Q-factor from the frequency-dependent PLDOS enhancement spectra, taking the metasurface with finite size 𝑁 = 5 as an example (see Figure R5). The first peak in the spectrum is identified as the fundamental mode of the BIC mode, while the higher-order peaks are attributed to collective hybridized modes formed through the coherent coupling of localized resonances across the unit cell array. These resonances increase in number and collapse into the BIC mode when the size of the array increases.\nTo extract the Q-factor, we performed a Fano line-shape fitting on the first resonant peak using the following equation in Eq.1,\n2\n(𝑞+\n2(𝑓−𝑓0)\n)\n𝑦 = 𝐴× ∆𝑓 +𝐵, (1)\n2\n1+(\n2(𝑓−𝑓0)\n)\n∆𝑓\nwhere 𝑓 is the resonant frequency of the peak, ∆𝑓 is full width at half maximum (FWHM), 0 and 𝑞 is Fano asymmetry parameter. In addition, A and B are scaling and offset constants, respectively. The quality factor is then calculated as: 𝑄 = 𝑓 /∆𝑓. This method provides an accurate characterization of the fundamental resonance associated with the BIC mode.\nFigure R5 Extracted Q-factor obtained by fitting the first resonant peak in the PLDOS enhancement spectrum using a Fano line profile. The example shown corresponds to a metasurface of finite size 5×5 .\nTo address the reviewer’s suggestion, we have now included a detailed description of the Q-factor extraction method in the Supplemental Information of the revised manuscript.\nFurthermore, we have included the experimentally measured evolution of the Q factor and PLDOS enhancement as a function of the number of dimers 𝑁×𝑁, as shown in Figure R6. Both quantities initially increase with array size, reflecting the growing coherence and collective response of the system. The PLDOS enhancement saturates (or grows very slowly ) as the metasurface becomes larger, in good agreement with simulation results. We note, however, that the experimentally retrieved Q values and PLDOS enhancement are lower than those obtained from simulations shown in Figures 3(f) and 3(g) (in the revised manuscript). This discrepancy arises primarily due to the limited frequency resolution of our measurement system (20 GHz), as well as the presence of material losses and substrate effects, which are not accounted for in the simulations. To better illustrate the ideal characteristics of the BIC, namely, the divergence of the Q factor and the saturation of the PLDOS with increasing metasurface size, we have chosen to present the simulated results in the main text. We sincerely appreciate the reviewer’s insightful comments, which have helped improve the clarity and completeness of our manuscript.\nFigure R6. The extracted Q-factors and PLDOS enhancement of the BIC mode from the measurement of the metasurfaces with finite sizes composed of different numbers of dimers.\n5. In Figure 5(a) and (b), the calculated Pyy/P0 exhibits intensity distribution below the BIC frequency, whereas the experimental results show prominent intensity above the BIC frequency. The authors should explain this discrepancy. Additionally, the authors should provide frequency-domain spectra at multiple height levels（the results of four height in Figure 5(a) and (b)）, Such data would more clearly illustrate how the PLDOS varies spatially and spectrally near the BIC frequency.\nResponse: Thank you for the question and the suggestion. The experimental results show a prominent intensity above the BIC frequency due to an incomplete subtraction of the direct interaction between the two microprobes in THz transient signal, which is from -10 to 6 ps. In the revised manuscript, it has now been corrected.\nRegarding the frequency spectra at multiple height levels corresponding to the THz transient, we have added these data in the revised manuscript, as shown in Figure R7. The following description has been added:\n‘To clearly show the decay of the BIC mode from the surface, the amplitude spectra are obtained by performing a Fourier transform of the THz transient signal, shown in Figure 5(a), after subtracting the contribution from direct dipole interaction. A clear peak at 0.395 THz which corresponds to the BIC mode is observed in the frequency spectra for both y- and x-polarized THz excitation and detection, respectively shown in Figure 5(b), and the intensity of the peak rapidly decreases with increasing probe height above the surface.’ (Page 12 in the revised manuscript )\nFigure R7. Out-of-plane confinement of the PLDOS enhancement for the BIC. (a) THz transients and (b) THz amplitude spectra measured by the near-field microprobes E and E . The measured points are (90 µm, 250 µm) y x for y-polarization and (180 µm, 210 µm) for x-polarization, corresponding to the position in Figure 2(b). (c-d) Measured PLDOS enhancements P /P and P /P as a function of frequency and height for y- and x- yy 0 xx 0 polarizations, respectively. (e-f) Calculated P /P and P /P at positions (90 µm, 250 µm) and (185 µm, 210 µm), yy 0 xx 0 respectively, as a function of frequency and height. (g) Decay of the PLDOS enhancement at the BIC frequency with h/λ , where h is the height above the metasurface and λ is the wavelength of the BIC mode, for (i) y- 0 0 polarized, and (ii) x-polarized measurements. The calculated (lines) and measured PLDOS (symbols) are normalized by the maximum value.\nReviewer #3 (Remarks to the Author):\nThe authors present a near-field THz-TDs that enables the direct spatial mapping of BICs in metasurfaces. The results provide the evidence of a large partial LDOS (PLDOS) enhancement associated with BICs. In my view, this work presents several intriguing aspects that are valuable to readers and colleagues in this field. However, there are several points that require clarification.\nResponse: We sincerely thank Reviewer #3 for the positive assessment of our work. In response to the reviewer’s specific comments below, we have carefully addressed each point by providing additional explanations, modifying the text for clarity, and including supplementary details where appropriate. We believe these revisions improve the clarity and rigor of the manuscript, and we hope they satisfactorily resolve the reviewer’s concerns.\n1. It should clearly explain the “nonlocal metasurfaces” to ensure that readers understand the nonlocal in this manuscript.\nResponse: To clarify the concept of “nonlocal metasurfaces” and highlight its relevance to our work, we have revised the manuscript to include a more detailed explanation. In resonant systems, high-Q resonances can be realized through long-range interactions, such as guided surface waves and surface lattice resonances, mechanisms that are foundational to nonlocal metasurfaces. Unlike local metasurfaces, where each unit cell responds independently to the incident field, nonlocal metasurfaces support optical modes that emerge from coherent interactions among many identical neighboring units. These collective interactions give rise to phenomena such as guided-mode resonances (GMRs) and bound states in the continuum (BICs), enabling high-Q flat optics.\nIn this work, the nonlocal metasurface consists of dimer arrays composed of two identically sized gold rods. The BIC mode arises from spatially extended interactions facilitated by surface lattice resonances, which indicates its nonlocal behavior. In contrast to idealized simulations, all fabricated metasurfaces used in the experiments have finite dimensions, which inherently reduce the achievable Q-factors compared to theoretical predictions. Moreover, the investigation of finite-size effects on Q-factors has predominantly been limited to theoretical studies, due to the challenges associated with traditional far-field measurement techniques. In this study, leveraging a near-field terahertz (THz) setup, we experimentally investigate the size-dependent Q-factor and PLDOS enhancement (corresponding to Figure 3 in the manuscript).\nWe have added the following description of nonlocal metasurfaces to the revised manuscript\n‘Unlike local metasurfaces, where each unit cell or meta-atom responds independently to the incident field, nonlocal metasurfaces support optical modes that arise from coherent interactions among many identical neighboring units. This nonlocality is evident in the behavior of the BIC mode in the metasurface system studied here.’ (Page 7-8 in the revised manuscript )\n2. The authors state that this is the first direct experimental evidence of LDOS enhancement by a BIC. However, in Ref. 24, an “optically inaccessible” BIC mode has been probed by electron energy-loss spectroscopy (EELS). They should compare this work to Ref. 24 in detail.\nResponse: We thank the reviewer for this helpful comment and for bringing Ref. 24 to our attention. We agree that Dong et al. (Ref. 24) provide an important experimental investigation of BIC modes using electron energy-loss spectroscopy (EELS) and cathodoluminescence (CL) in a scanning transmission electron microscope (STEM). Below, we provide a detailed comparison of their work and ours to clarify the complementary nature of these two approaches.\nFirstly, both studies aim to assess the optical properties of symmetry-protected BICs, which are typically inaccessible through conventional far-field optical measurements. Dong et al. use a high-energy electron beam to locally excite BIC modes in silicon nanoantenna arrays, probing the system’s photonic states via EELS and CL. Their approach provides nanoscale spatial resolution and is particularly powerful for mapping BIC modes in the visible and near-infrared spectral range. In contrast, our work uses a near-field THz time-domain spectroscopy (THz-TDS) setup with sub-wavelength photoconductive microprobes to directly excite and detect the near-field optical response under THz electromagnetic excitation. This allows us to map the partial local density of optical states (PLDOS) in real space and observe its spectral and spatial evolution as a function of metasurface size. Our measurements provide direct evidence of PLDOS enhancement associated with BICs, under electromagnetic (rather than electron-beam) excitation conditions.\nSecondly, Dong et al. demonstrate that BICs can be accessed in the near-field using EELS, providing experimental evidence of their spatially confined character. The measured signal reflects the electron-induced excitation spectrum, rather than the field distribution that would result from direct electromagnetic excitation. In our work, the measured PLDOS is obtained under relevant condition to light–matter interaction in photonic applications, such as spontaneous emission and dipole coupling, and our data reveal not only the spatial localization of the BIC but also the saturation behavior of the PLDOS despite a diverging Q-factor.\nIn the revised manuscript, we have included an explicit reference to Dong et al. in the context of prior work, acknowledging its contribution and clarifying its distinctness as below:\n‘In related work, Dong et al. used electron energy-loss spectroscopy (EELS) to probe the near-field response of optically inaccessible BIC modes in silicon nanoantenna arrays, revealing strong spatial confinement of the associated electromagnetic energy with nanometer resolution.’ (Page 3 in the revised manuscript )\n3. The leakage of BIC states is generally observed by the quasi-BIC patterns. The author should more clearly present the influence of the probe on the BIC/quasi-BIC mode.\nResponse: We appreciate the Reviewer’s valuable comment. We agree that the influence of the probe on the BIC/quasi-BIC mode needs a clearer explanation. This measurement technique accurately determines the change of the PLDOS due to the photonic surroundings of the source. The probes do not alter this PLDOS significantly, which has been explored in one of our previous works (Huurne et al. ACS Photonics 10, 2980 (2023)). Figure R1 reproduces a figure of this manuscript. Figure R1(a) shows a map of the imaginary component of the electric field amplitude as a function of frequency and distance to a gold mirror relative to the wavelength. The clear oscillations of the PLDOS are visible as a function of the distance between the probes and the interface. The oscillatory behavior of the PLDOS is characteristic of the self-action of the emitted radiation by the point dipole on itself due to the reflection at the interface. In addition, Figure R1(b) compares the measured PLDOS of gold, InSb and quartz interfaces with the respective theoretical calculation for a point dipole as a function of height at 0.76 THz. Note that there are no fitting parameters in these calculations. The excellent agreement between measurements and theory demonstrate that the probes do not alter significantly the PLDOS.\nFigure R1. (a) Imaginary component of the normalized Fourier transform of the electric field amplitude as a function of frequency and distance to the mirror relative to the wavelength. Comparison between the measured PLDOS enhancement and analytical result at 0.76 THz as a function of distance from an interface of (b) gold, (c) InSb, and (d) quartz.\nTo clarify this point, we have added the following sentences in the revised manuscript:\n‘Based on this unique setup, we have previously investigated the PLDOS enhancement above gold, InSb, and quartz interfaces, achieving excellent agreement with theoretical calculations for a point dipole as a function of height and demonstrating that the microprobes do not alter significantly the PLDOS [20]. (Page 6 in the revised manuscript)\n4. In Fig. 4, the authors severally present calculated and measured spatial maps of PLDOS enhancement. Are these maps all normalized to the same maximum?\nThank you for raising this important point, and we apologize for any confusion caused. The normalization procedures for the simulated and experimental PLDOS enhancement maps are indeed different. Specifically, the calculated PLDOS enhancements are normalized to the corresponding values in vacuum at each frequency, whereas the experimental results are normalized to the PLDOS measured on the bare quartz substrate at the same frequency.\nWe have clarified this normalization approach in the revised manuscript by adding the following description:\n“The calculated PLDOS enhancements are normalized to the corresponding values in vacuum at each frequency, whereas the experimental results are normalized to the PLDOS measured on the bare quartz substrate at the same frequency.” (page 10 in the revised manuscript)\n5. The authors specifically explore the enhancement of PLDOS depending on different numbers of dimers N. whether the other parameters, for example, the thickness of gold film, will affect the enhancement?\nResponse: We thank the reviewer for this valuable question. In our study, we focused on the influence of metasurface array size (i.e., the number of dimers) on the PLDOS enhancement, while keeping structural parameters such as rod dimensions, spacing, and gold thickness fixed. The thickness of the gold film is an important factor that can influence the PLDOS through several mechanisms. If the film is too thin (i.e., comparable to or less than the skin depth at THz frequencies), it may not fully support confined surface currents, leading to increased ohmic losses and reduced reflectivity. This increased material absorption can degrade the Q-factor of the resonance, thereby reducing the PLDOS enhancement. The PLDOS is sensitive to how tightly the electromagnetic fields are confined within the resonator. Thicker gold rods generally improve confinement. Conversely, thinner rods may result in weaker confinement and broader field profiles, increasing the effective mode volume and reducing the PLDOS enhancement.\nThe efficiency with which individual resonators couple to each other depends on their geometrical cross-section, including thickness. Insufficient thickness may weaken the scattering strength of each dimer, thereby impairing the formation of coherent collective modes such as the symmetry-protected BIC.\nFollowing on the comment of Reviewer #3, we have added the following:\n‘While our study focuses on the influence of array size on PLDOS enhancement with fixed structural parameters, we note that other parameters, such as the gold film thickness, can also significantly affect the PLDOS. This occurs through its impact on material losses, field confinement, and collective mode formation’ (Page 13 in the revised manuscript)"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have thoroughly addressed the reviewers' concerns. In my opinion, the manuscript can be published in Nature Communications.", "reviewer_2": "The authors have provided detailed responses to the various questions raised by the reviewers. I believe this article demonstrates the enhancement of the photonic density of states by dimers from a near-field perspective, but the presentation must be more rigorous, as I do not think this work has proven that the various hybrid resonances are BIC resonances. The excitation conditions for BIC are stringent; to ensure it is above the light line, theoretical BIC states require normal incidence (corresponding to the Γ point in the band structure). The excitation method used in this work makes it difficult to believe that only pure BIC states were excited. Furthermore, from a detection perspective, the symmetry-protected BIC states of dimers are highly antisymmetric, and I suspect that the proximity of the probe could directly disrupt the BIC state, causing it to leak into quasi-BIC states or other states. The experimental validation approach adopted in this paper also fails to prove that the detected resonances are BIC states: First, Fig. 3e clearly shows varying brightness among different rows of dimers, resembling propagating traveling-wave modes along the surface. How can the authors be certain that the excitation is not due to discrete states at the M or X points below the light line? Second, compared to Fig. 4b, the measured field distribution in Fig. 4e lacks the expected symmetry in both the Y and X directions. Crucially, can the authors provide phase information to demonstrate the antisymmetric relationship between the two bars within a dimer? Third, regarding the disappearance of the resonance signal when the probe is lifted to a certain height, I do not believe this proves the detection of a state confined to the surface. It seems more likely due to weak point-source excitation—that is, the signal excited by the point-source probe close to the metasurface would naturally attenuate as the probe moves away. This is probably a natural consequence of increasing the distance between the two probes and does not prove the detection of a confined state.\nIn summary, this work validates the enhancement of the photonic density of states by resonant states from a near-field perspective, but it should not emphasize that the enhancement originates from BIC states. Statements in the abstract such as, \"With a unique terahertz near-field microscope, we have mapped the spatial distribution of BICs in nonlocal metasurfaces,\" are particularly unfounded. If the authors can revise the abstract, conclusion, and full text to avoid giving readers the impression that they detected pure BIC states, and cautiously describe that the cause of the enhancement in the photonic density of states should be attributed to quasi-BIC states rather than pure BIC states (which may be included but are unlikely to be the only states), then I believe the editor may reconsider accepting this article.", "reviewer_3": "The manuscript has been revised according to suggestions, and we recommend it for publication."}, "author_response": "Response letter\nReviewer #1 (Remarks to the Author): The authors have thoroughly addressed the reviewers' concerns. In my opinion, the manuscript can be published in Nature Communications.\nReviewer #3 (Remarks to the Author): The manuscript has been revised according to suggestions, and we recommend it for publication.\nResponse: We sincerely thank Reviewers #1 and #3 for their positive evaluation and for recommending the acceptance of our manuscript. We greatly appreciate the time and effort that the reviewers have dedicated to reviewing our work and for the constructive suggestions in earlier rounds, which have helped us improve the quality and clarity of the manuscript.\nReviewer #2(Remarks to the Author):\nThe authors have provided detailed responses to the various questions raised by the reviewers.\nI believe this article demonstrates the enhancement of the photonic density of states by dimers from a near-field perspective, but the presentation must be more rigorous, as I do not think this work has proven that the various hybrid resonances are BIC resonances. The excitation conditions for BIC are stringent; to ensure it is above the light line, theoretical BIC states require normal incidence (corresponding to the Γ point in the band structure). The excitation method used in this work makes it difficult to believe that only pure BIC states were excited. Furthermore, from a detection perspective, the symmetry-protected BIC states of dimers are highly antisymmetric, and I suspect that the proximity of the probe could directly disrupt the BIC state, causing it to leak into quasi-BIC states or other states. The experimental validation approach adopted in this paper also fails to prove that the detected resonances are BIC states:\nResponse: We thank the reviewer for raising the point regarding the rigorous identification of BIC resonances. We agree with the reviewer that the distinction between symmetry-protected BICs and quasi-BICs is subtle.\nThe metasurface under study supports a symmetry-protected BIC (commonly associated with the Γ point in the Bloch band structure). Under plane-wave excitation at normal incidence, such modes remain fully uncoupled from the far field due to symmetry mismatch and are therefore invisible in the transmission spectrum, as shown in Fig. 1. Our near-field experimental approach circumvents this restriction by employing two microprobes that act as a localized dipolar source and detector. The localized probes can excite and detect the (quasi-)BIC mode within the unit cell, bypassing the symmetry mismatch that prevents the far-field excitation. As pointed out by the reviewer, the emitter and detector probes locally break the symmetry of the metasurface, and allows the coupling to the BIC. We should stress that this effect is very local due to the small size of the microprobes compared to the wavelength. Therefore, this locally broken symmetry, needed to couple in and out the metasurface, is much less pronounced than the broken symmetry of metasurfaces with asymmetric unit cells defining quasi-BICs, as commonly discussed in the literature. In this last case, the symmetry is broken in each unit cell (instead of only at the excitation and detection spots). Therefore, it is impossible to excite pure BICs in its most strict definition. However, our large metasurface, formed by the symmetric array of dimers of rods with equal length, and the measuring technique with local probes for excitation and detection make this system the closest approximation to a true BIC investigated so far.\nOne important remark is that pure symmetry-protected BICs are only possible in infinite systems, which are obviously impossible to attain in any experiment. Boundaries in finite systems break the translation symmetry and introduce a loss channel. Therefore, the the small-sized arrays that we also investigate deviate from the idealized case.\nWe have clarified these points in the revised manuscript, stressing the fact that true (theoretical) BICs in their strict definition are experimentally impossible to realize due to the inherent finite size of samples and to measure due to the unavoidable need of locally breaking the symmetry of the lattice. Following the criticism of the reviewer and to avoid any possible misinterpretation of our results, we have labelled the measured modes as quasi-BICs.\n1. First, Fig. 3e clearly shows varying brightness among different rows of dimers, resembling propagating traveling-wave modes along the surface. How can the authors be certain that the excitation is not due to discrete states at the M or X points below the light line?\nResponse: We thank the reviewer for this thoughtful question regarding the interpretation of the spatial patterns in Fig. 3(e). By definition, a traveling wave propagates through space and time, with points continuously advancing in phase, whereas a standing wave, typically formed by the superposition of counter-propagating waves, oscillates in place with fixed nodes and antinodes. The brightness variations across different rows in Fig. 3(e) resemble standing-wave–like patterns, which arise from the finite-size discretization of the Γ-point BIC mode. In a finite metasurface, the symmetry-protected BIC of the infinite array hybridizes into a set of discrete collective modes confined within the array boundaries. These hybridized modes manifest as intensity modulations across rows or columns, but unlike traveling modes, they correspond to standing-wave field distributions determined by the finite array size. This interpretation is consistent with the sharp resonances in Fig. 3(d), which split into multiple discrete peaks as the array size increases.\nRegarding the possibility of discrete states at the M or X points below the light line, we emphasize three points. Firstly, M/X states exhibit characteristic π phase alternations between neighboring unit cells (checkerboard or alternating-column patterns), which are not observed in our measurements. Instead, the measured distributions repeat in phase across rows and columns, consistent with Γ-point modes. Secondly, the observed out-of-plane confinement length (~λ₀/20) is much shorter than the ~0.2λ₀ expected for X/M states, further ruling out their contribution. Finally, M/X points typically occur at different frequencies due to band dispersion. In our case, simulations confirm that M/X modes are not supported by the system because of the substrate.\n2. Second, compared to Fig. 4b, the measured field distribution in Fig. 4e lacks the expected symmetry in both the Y and X directions. Crucially, can the authors provide phase information to demonstrate the antisymmetric relationship between the two bars within a dimer?\nResponse: We thank the reviewer for this insightful observation. The expected symmetry in both the Y and X directions that the reviewer refers to is not visible in the measurements because emitter and receiver need to be displaced by a small distance of about 10 mu to avoid cross-talking due to the laser light separately reflected on the emitter photoconductive switch reaching the detector. This small displacement causes the asymmetry as described on page 10 of the revised manuscript.\nFigures 4b and e display the measurements of the PLDOS. Therefore, in these measurements, it is not expected to see the antisymmetric relationship between the two bars within a dimer that the reviewer mentions. The PLDOS in enhanced on all the edges of the rods. The antisymmetric phase relation is expected only for the surface currents and fields’ amplitudes. This antisymmetric phase relation for the fields has been demonstrated earlier by some of us in the manuscript “Unveiling the Symmetry Protection of Bound States in the Continuum with Terahertz Near-Field Imaging, ACS Photonics 2021, 8, 3010−3016”. Figure 6 of this manuscript is reproduced below, showing the measurements and calculations for the z-field component (out-of-plane field). Please, note that the emitter in these measurements is located at position (0,0), while the receiver is scanned over the array. The noise along the direction X~0 is due to the cross-talk between emitter and receiver due to the laser reflected on the emitter reaching the detector (the cross talk mentioned above).\nWe have confirmed these measurements in our experimental configuration, as illustrated in the newly added Section S6 of the Supplementary information (Near-field measurements of the phase). In this section, we display spatial coherence measurements over one unit cell by showing the imaginary component of the y- and z-field components. In these near-field measurements, the antisymmetric fields between the two rods in the dimer are clearly visible. We believe that this addition to the supplementary information demonstrates the nature of the quasi-BIC and improves significantly the manuscript. We are thank full to the reviewer for motivating these improvements.\n3. Third, regarding the disappearance of the resonance signal when the probe is lifted to a certain height, I do not believe this proves the detection of a state confined to the surface. It seems more likely due to weak point-source excitation—that is, the signal excited by the point-source probe close to the metasurface would naturally attenuate as the probe moves away. This is probably a natural consequence of increasing the distance between the two probes and does not prove the detection of a confined state.\nResponse: We investigate the partial local density of optical states (PLDOS) resulting from resonant photonic environments, specifically the symmetry-protected quasi-BIC mode, which is invisible from far-field measurement and spatially confined near the metasurface surface. By definition, the PLDOS is obtained from the field at the position of the source, and consequently, the emitter and detector are very close to each other at a fixed distance. It is important to realized that in our measurements, the two microprobes move together away from the surface.\nTo further confirm that the resonance corresponds to a surface-confined mode rather than trivial attenuation of a point-source signal, we have performed an additional height-dependent measurement in which the emitter remained fixed relative to the metasurface, while the detector probe was translated away from the sample in 2-µm steps, as illustrated in Figure. R2(a). At each height, we recorded the THz transient. The 2D map in Figure R2(b) shows a clear field enhancement at the BIC frequency (~0.395 THz), which rapidly disappears within ~25 µm above the surface. The extracted resonance amplitude versus height is plotted in Fig. R2(c). The measured data (symbols) are well described by an exponential decay function, 𝑦 = 𝑦 + 0 𝐴𝑒𝑥𝑝ቀ ℎ−ℎ0ቁ, where 𝜏 is the decay length. The fitted decay constant corresponds to 𝜏 ≈ 𝜏 0.032𝜆 (≈ 20 µm for 𝜆 ≈ 760 µm). This sub-wavelength confinement length is consistent 0 0 with the evanescent nature of a symmetry-protected BIC-derived mode, and is in excellent agreement with the rapid PLDOS decay observed in our main measurements. By contrast, a simple reduction of point-source coupling with probe height would lead to a broadband attenuation rather than the sharp, narrow-frequency exponential decay we observe. Together, these results confirm that the disappearing signal with increasing height in Fig. 5 in page 12 is a direct signature of the spatial confinement of the BIC-derived mode to the metasurface.\nFigure R2. The out-of-plane confinement of BIC mode to the surface. (a) Illustration of the measurement setup; (b) 2D colorplot of the measured THz field intensity as functions of frequencies and height between the surface and the detected probe. (c) The normalized THz intensity as a function of height at the BIC mode. The symbol are the measured data, and the line is the fitted result with an exponential decay function.\nThese new experimental results on BIC mode’s high confinement to the surface have been added in Section S8 of Supplementary Information (Out-of-phase confinement of the quasi-BIC mode to the surface).\nIn summary, this work validates the enhancement of the photonic density of states by resonant states from a near-field perspective, but it should not emphasize that the enhancement originates from BIC states. Statements in the abstract such as, \"With a unique terahertz near-field microscope, we have mapped the spatial distribution of BICs in nonlocal metasurfaces,\" are particularly unfounded. If the authors can revise the abstract, conclusion, and full text to avoid giving readers the impression that they detected pure BIC states, and cautiously describe that the cause of the enhancement in the photonic density of states should be attributed to quasi-BIC states rather than pure BIC states (which may be included but are unlikely to be the only states), then I believe the editor may reconsider accepting this article.\nResponse: Following the reviewer's recommendation, we have revised the abstract and conclusions to emphasize that the investigated metasurfaces will only support a true BIC in their infinite (idealized) counterparts. While in our measured samples, due to the inherent finite size and local broken symmetries introduced by microprobes, we do not have pure BICs. In the revised manuscript, We have corrected it as quasi-BIC in the maintext, which is marked in red.\nFurthermore, the revised abstract and conclusion is shown as below.\nNew abstract:\nBound States in the Continuum (BICs) are exotic optical modes that remain decoupled from free-space radiation. Symmetry-protected BICs in infinite metasurfaces provide unique opportunities to tailor light–matter interactions in extended photonic structures. One of the most important properties of these structures is the partial local density of optical states (PLDOS), which describes the number of electromagnetic modes available for a photon to occupy at a specific position, frequency, and polarization. It quantifies how the surrounding photonic environment modifies the ability of light to couple with matter at a given location. Here, we employ a terahertz near-field microscope with dual local probes to directly excite and detect quasi-BICs in finite metasurfaces with inversion symmetry, which support a symmetry-protected BIC in their infinite counterpart. Our results provide the first direct evidence of a large PLDOS enhancement associated with quasi-BICs in finite metasurfaces. As the metasurface size increases, the quasi-BIC evolves into the BIC, and the quality factor (Q) diverges while the PLDOS saturates at a finite value. This result reflects an upper limit for enhanced light-matter interaction by BICs. Our findings pave the way for next-generation on-chip metasurfaces with maximum light-matter interaction strengths.\nNew conclusions:\nIn this work, we have used two THz microprobes positioned at a subwavelength distance from each other to locally excite and probe the complex electric field distribution of modes in metasurfaces with inversion symmetry formed by arrays of Au dimer rods. The infinite size counterpart of these metasurfaces support symmetry protected BICs. The near-field probes allow the direct measurements of the partial density of optical states of quasi-BICs associated with the finite metasurfaces, which was found to be significantly enhanced at the edges of the rods, extending over the unit cell due to the collective resonant field enhancement. Numerical calculations show the saturation of the PLDOS with the size of the array at moderate sizes of a few unit cells, despite the increase in the Q-factor for larger arrays. The increase in Q-factor in symmetry-protected quasi-BICs is unavoidably linked to the increase of the mode volume due to the increase of the array’s size and the nonlocal character of these modes, which leads to the saturation of the ratio Q/V that defines the PLDOS enhancement. Furthermore, we theoretically and experimentally investigated the spatial confinement to the surface of the PLDOS enhancement at the quasi-BIC frequency, revealing a strong confinement accompanied by a large enhancement.\nOwing to the scale invariance of Maxwell’s equations, the design principles and insights developed here for terahertz metasurfaces can be naturally extended to the infrared and visible regimes. In particular, the ability to engineer (quasi-) BIC modes and enhance the PLDOS opens up promising avenues for controlling light–matter interactions in nanophotonic systems. This has important implications for photonic and quantum technologies, such as the enhancement of the emission efficiency of single-photon sources, the tailoring of the spontaneous emission rates in cavity quantum electrodynamics, or the improvement of nonlinear optical responses in integrated photonic platforms."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The author's response has largely addressed my core concern in the review, namely that what was detected is a quasi-BIC rather than a pure BIC. Given the positive implications of this work for the study of near-field BIC metasurfaces, I recommend acceptance."}, "author_response": ""}]

Decision letter