Ultrathin 3R-MoS2 metasurfaces with atomically precise edges for efficient nonlinear nanophotonics

Abstract

Dielectric metasurfaces that combine high-index materials with optical nonlinearities are recognized for their potential in quantum and classical nanophotonic applications. However, the fabrication of high-quality metasurfaces poses material-dependent challenges, as their designs are often susceptible to disorder, defects, and scattering losses, mostly occurring at the edges of nanostructured features. Additionally, the choice of the material platforms featuring second-order optical nonlinearities χ⁽²⁾, is limited to broken-inversion symmetry crystals such as GaAs, GaP, LiNbO₃, and various bulk van der Waals materials. Here, we use a combination of top-down lithography and anisotropic wet etching of a specially stacked van der Waals crystal – 3R-MoS₂, which exhibits both a high refractive index and exceptional χ⁽²⁾ nonlinearity, to produce ultrathin (~20–25 nm) metasurfaces with atomically sharp edges, where the etching breaks the in-plane symmetry of the meta-atoms. The broken symmetry manifests as a quasi-boundstate-in-the-continuum, enabling the enhancement of second-harmonic generation of three orders of magnitude at specific wavelengths.

Dielectric metasurfaces with strong nonlinearities hold promise for integrated nanophotonics and quantum technologies, yet their fabrication is hindered by material-dependent challenges and limited material choices with second-order optical nonlinearities. Here, the authors employ top-down lithography and anisotropic wet etching to fabricate atomically precise, ultrathin metasurfaces highly efficient for second-harmonic generation.

Introduction

Transition metal dichalcogenides (TMDs) form a class of 2D materials widely recognized for their remarkable optical and electronic characteristics¹. Due to their unique crystalline structure and direct band gaps in the visible and near-infrared spectral range, 2H-TMD monolayers, with a general formula MX₂ (M = Mo, W and X = S, Se or Te), exhibit broken inversion symmetry, making them useful for applications in nonlinear nanophotonics and optoelectronics^{2, 3, 4–5}. Large values of the second-order nonlinear susceptibility ∣χ⁽²⁾∣ are essential for processes, like second-harmonic generation (SHG)⁶, sum- and difference-frequency generation (SFG, DFG)^7,8, optical parametric amplification (OPA)⁹, spontaneous parametric down-conversion (SPDC)^10,11, electro-optical effect (EO)¹², and others. Indeed, non-centrosymmetric TMD monolayers are characterized by impressive ∣χ⁽²⁾∣ values, reaching, for example, for MoS₂ ~10⁻⁷ m⋅V⁻¹¹³ (when expressed in bulk-equivalent units). Despite large ∣χ⁽²⁾∣ values, TMD monolayers face a challenge in achieving optimal nonlinear conversion efficiencies for SHG and SPDC, η^10,14, due to nanometric thicknesses and, therefore, wave propagation length (η_SHG = I_2ω/I_ω ∝ ∣χ⁽²⁾∣²L²¹⁵, where L is the propagation length). This is in contrast to conventional phase-matched bulk crystals that are able to frequency-double light with efficiencies exceeding 50%¹⁶.

A popular strategy to increase the nonlinear conversion efficiency of TMDs involves integrating the monolayer onto a resonant metasurface spectrally aligned with the active material’s peak ∣χ⁽²⁾∣ response^{17, 18, 19, 20–21}. Although this approach can enhance the nonlinear response, it does not directly address the volume constraints of TMD materials. Furthermore, it introduces additional complexities due to the involvement of semiconductor interfaces and strain effects during the exfoliation and integration of monolayers with resonant metasurfaces. These factors can result in spectral shifts of exciton resonances^22,23 and damping of nonlinear responses^24,25, negatively impacting the coupled system’s performance.

Hence, innovative approaches are required to unlock the full potential of various nanostructured materials for integrated nanophotonics and quantum technologies. In this regard, high-quality factor (high-Q) resonances in traditional dielectric nonlinear materials such as LiNbO₃^26,27, BaTiO₃²⁸, GaAs²⁹, GaP³⁰ with considerably smaller ∣χ⁽²⁾∣ values (~10⁻¹¹−10⁻¹⁰ m⋅V⁻¹) in comparison to monolayer TMDs are especially promising^{31, 32–33}. Various designs including metasurfaces^34,35, waveguides³⁶, and whispering gallery modes (WGM)³⁷ exhibit a dramatic efficiency increase of optical nonlinear processes including SHG and SPDC³⁸. The goal of this nonlinear all-dielectric nanophotonic concept is to identify materials with a high refractive index (n) that are robust to high-intensity excitation, have low optical losses (k), and high ∣χ⁽²⁾∣ nonlinearity to efficiently confine electromagnetic fields within resonant nanostructures and generate a strong nonlinear response.

A powerful contemporary concept to achieve optimal field confinement in nanophotonic structures is optical bound-state-in-the-continuum (BIC)^39,40, which benefits from high-index dielectric materials⁴¹. BICs feature electromagnetic field decoupling from the environment and, as a result, theoretically, an infinite radiative quality factor, which, however, in real structures is limited by fabrication intolerances and material losses^42,43. Such high-Q BIC resonances are challenging to excite from the far field due to optical reciprocity – low optical leakage leads to weak optical excitation. By introducing a slight asymmetry in the system, one can efficiently excite these high-Q resonances, at the expense of reducing their Q-factor⁴⁴. Such a concept, often referred to in the literature as quasi-BIC (q-BIC), has great potential in a variety of nonlinear processes. It has been recently demonstrated that the SHG conversion efficiencies in structures featuring q-BICs can reach  ~0.1%^{45, 46–47}, considerable SPDC^{48, 49–50}, and even high-harmonics generation^51,52.

Multilayer TMD materials are widely recognized for their high index (n > 4) and extremely anisotropic (Δn > 1.5) optical properties^53,54, which makes them attractive for the realization of BIC concepts. Specifically, their use encompasses enhanced light-matter coupling^55,56, leveraging (q-)BICs to reach the strong coupling regime^57,58, Mie modes to enhance the nonlinear signal originating from the broken inversion symmetry at the surface^59,60, and lasing in WGM disks⁶¹. However, conventional 2H phase bulk van der Waals crystals are centrosymmetric, resulting in a vanishing nonlinear signal^59,60. In contrast, alternative TMDs, such as GaSe, NbOCl₂, 3R-MoS₂, and 3R-WS₂, which possess substantial ∣χ⁽²⁾∣ nonlinearities even in bulk crystals (~10⁻¹⁰−10⁻⁹ m⋅V⁻¹), have recently attracted research attention, particularly in terms of nonlinear optical and piezoelectric applications^{6,11,62, 63, 64, 65, 66, 67, 68, 69, 70, 71–72}. However, to date, only a few high-index metasurfaces fabricated using these materials have been reported⁷³.

In this work, we demonstrate a bulk-χ⁽²⁾-active all-TMD metasurface supporting q-BIC resonances. Our high-index nonlinear 3R-MoS₂ metasurface with a deeply subwavelength 20–25 nm thickness is able to significantly enhance the SHG signal – more than three orders of magnitude compared to that of an optimized host flake. We effectively control the optical losses of the BIC system through a wet etching process applied to convert a circular hole metasurface into the one composed of perfect equilateral triangles with atomically sharp edges⁷⁴. This precise etching technique allows us to break the C₂ and C₄ in-plane symmetries, facilitating the excitation of the q-BIC resonances by far-field illumination. Furthermore, the ultrasharp nature of meta-atoms’ edges enable the experimental realization of high-Q optical resonances⁷⁵, as schematically illustrated in Fig. 1. Together with the recently developed wafer-scale fabrication methods of the 3R phase van der Waals materials⁷⁰, this all-TMD nonlinear high-Q concept holds promise to enhance the efficiency of a wide range of second-order nonlinear processes in deeply subwavelength nanophotonic settings.

[See PDF for image]

Fig. 1

Concept of broken-symmetry 3R-MoS₂ metasurface.

Schematic of the wet etching process which forms atomically sharp equilateral triangular holes from the pre-etched circular holes of smaller size obtained via reactive ion etching of 3R-MoS₂ multilayer flake. Wet-etching provides atomically sharp edges, clean surface, and broken C₂ and C₄ symmetry, enabling the realization of quasi-bound-state-in-the-continuum (q-BIC) concept. The square insets are SEM images of the typical circular and triangular hole unit cells. Note that since wet etching consumes part of the material, it is necessary to start with a smaller initial circular hole diameter to achieve the desired size of the final triangular hole. Therefore, to achieve a comparable filling factor (ff) between circular and triangular metasurfaces, the initial circular hole size in the latter must be smaller in comparison to that of the former (see the dotted line representing the size of the original circular hole, as shown in the scanning electron microscopy image of the triangular nanohole). Scale bars of the SEM insets are 200 nm.

Results and discussion

3R-MoS₂ flake thickness optimization

We begin by optimizing the 3R-MoS₂ material platform through analysis of the SHG signal as a function of both pump wavelength and the thickness of unpatterned flakes. This approach allows us to gain insight into how optical resonances within a metasurface can amplify the inherently efficient nonlinear response observed in bare flakes. Although similar studies have been previously reported, they focused on a specific or limited range of pump wavelengths^{63, 64–65}. Our work aims to expand this understanding by exploring a wider range of wavelengths and preparing a solid foundation for following metasurface experiments in reflection geometry.

Figure 2 shows an experimental and theoretical study of the total intensity of the SHG as a function of a pump wavelength (λ_pump in 800–1040 nm range or photon energy ω_pump in 1.2–1.55 eV range) and the height of the bare 3R-MoS₂ flake exfoliated onto a glass substrate. The schematic of such a study is illustrated in Fig. 2a, where a tunable femtosecond laser is focused through an objective on the 3R-MoS₂ flake surface from the glass substrate side. The SHG signal is then collected through the same lens and quantified by an avalanche photodiode (APD), provided that SHG spectra are free of other signals⁶⁸. As shown in Fig. 2b, (i) an optical image of the flake is correlated with (ii) SHG mapping of the same flake at different pump wavelengths (example of a single wavelength SHG mapping at a single pump wavelength; SHG spectra, as well as SHG signal vs. pump power are presented in Supplementary Information (SI), Supplementary Fig. S1). Finally, (iii, iv) atomic force microscopy (AFM) mapping scans and profiles are measured to assess the relevant flake thicknesses. This approach allows us to collect SHG spectra for several flake thicknesses simultaneously.

[See PDF for image]

Fig. 2

Flake thickness optimization for second harmonic generation (SHG).

a Sketch of the SHG measurement of the 3R-MoS₂ flakes exfoliated on glass substrate. The height of each flake was measured with atomic force microscopy (AFM). b Standard study steps of the flake analysis: (i) Optical image of the flake of interest. Scale bar is 10 μm; (ii) SHG mapping of the intensity at a given pump wavelength; (iii) and (iv) AFM analysis of the flake height. The white dashed line in (i) corresponds to the AFM area of analysis in (iv). c Experimental SHG intensity of 3R-MoS₂ flakes for various pump wavelengths and flake thicknesses. The flake thicknesses, h, are measured by AFM. Dashed lines correspond to cross-sections in (d). d SHG intensity cross-sections at given pump wavelengths. Markers represent experimental points. Solid lines connect neighboring points for visual perception. The error bars in y direction are associated with the APD sensitivity. The error bars in x direction are associated with the AFM precision. Both of them are less than the marker size. e Theoretically calculated SHG intensity of the 3R-MoS₂ flakes using the transfer matrix method. Dashed lines correspond to cross-sections shown in (f). f SHG intensity from flake thickness at given pump wavelengths.

The experimental result of the SHG intensity of the 3R-MoS₂ flakes for different pump wavelengths, λ_pump, and flake thicknesses, h, is shown in Fig. 2c. Clear bands with high-intensity SHG are observed. They disperse toward increased thickness with decreasing pump energy, indicating resonant Fabry-Pérot (FP) behavior. The bright yellow parts of the map denote the maximum local SHG signal, which is determined by both the resonant FP conditions and the spectral dependence of the experimental χ⁽²⁾ tensor. Here, the data show a narrow local SHG maximum at small h (~20 nm), whose amplitude is slightly smaller than that at larger thicknesses, which agrees with previous data for reflection SHG spectroscopy of 3R-MoS₂ flakes^{63, 64–65}. This peak is attributed to the product of the transmission through the flake and the phase mismatch factor according to the following equation¹⁵:

1

$\frac{I_{2\omega }}{I_{\omega }}\propto \frac{T\left(\omega ,h\right)\mid {\chi }^{\left(2\right)}\left(\omega \right){\mid }^2{I}_{\omega }}{n_{2\omega }{n}_{\omega }^2}{h}^2{\mathrm{sinc}}^2\left(\frac{\mathrm{\Delta }kh}{2}\right)$where T(ω, h) is the transmission through a 3R-MoS₂ slab of thickness h at frequency ω, $\mathrm{\Delta }k=k_{2\omega }-2k_{\omega }=4\pi {\lambda }_{\mathrm{pump}}^{-1}\left(n_{2\omega }-n_{\omega }\right)$ – is the phase mismatch between the second-harmonic and pump beams, λ_pump is the pump wavelength, and n_2ω and n_ω are the refractive indices at the second-harmonic and pump wavelengths, respectively. Transmission peaks at h → 0, while the phase factor peaks at h_coh = π(Δk)⁻¹ > 0, with their overlap resulting in a local maximum at h ~ 20–25 nm, depending on the wavelength, which coincides with the first optimal 3R-MoS₂ flake thickness (Fig. 2c–f).

The relative amplitudes and widths of the SHG signal are further compared in Fig. 2d for selected pump wavelengths (marked in Fig. 2c by dashed lines) λ_pump: 800 nm (blue) and 1000 nm (black) – as the limits of the spectral range of the pump laser, and 910 nm (red) as the most resonant wavelength. This head-to-head comparison illustrates the sensitivity of the FP resonances to the exact value of h, providing up to three orders of magnitude enhancement in SHG intensity compared to off-resonant flake thicknesses.

The theoretical analysis of the SHG intensity was carried out using the transfer-matrix method (TMM) of a glass/3R-MoS₂ interface utilizing the χ⁽²⁾ tensor calculated in ref. ⁶⁸ and refractive index from ellipsometric measurements⁵⁴ (for more detailed information see Methods). Figure 2e presents calculated SHG for different heights of the flake (up to 400 nm) in the same energy range as in experiments. The dispersive FP resonances visible in the calculations match very well the experimentally observed dispersion. However, local SHG maxima within each FP band are slightly offset to the red versus the experiments, showing maxima at around 1.27 eV (~975 nm), where the theoretical χ⁽²⁾ tensor peaks⁶⁸. The dashed lines for a cross-section comparison are selected at different pump photon energies due to a slight spectral mismatch between the experimental and theoretical χ⁽²⁾ tensors. However, the behavior seen in Fig. 2f is similar to that observed experimentally, with a  ~3-order of magnitude variability in SHG for a given pump wavelength.

A comparison of the experimental and calculated data underlines that FP interference is the main factor responsible for the observed dispersion of the SHG signal, while the χ⁽²⁾ tensor modifies the spectral intensity along each FP band. Furthermore, the calculated SHG signal shows saturation with increasing thickness of the 3R-MoS₂ material for higher-order FP modes. This is most likely caused by the absorption at the second-harmonic wavelengths, where λ_SHG = λ_pump/2. Indeed, the imaginary part of MoS₂’s dielectric function in the 400–500 nm range exceeds 10⁵⁴, preventing light emission from deeper layers of the TMD. While a significant amplitude of the second-order polarization exists throughout the bulk of the flake, only part of it near the surface can radiate efficiently into the far field, while the remainder is strongly attenuated. Hence, we conclude that  ~20–25-nm-thick 3R-MoS₂ flakes are the most promising in terms of SHG emission per amount of material used in the studied spectral range. The choice of an ultrathin optimal thickness is further supported by several factors, including ease of nanofabrication, the absence of stacking faults (as discussed later), reduced optical losses (especially at the second-harmonic wavelength), and the lack of phase-matching constraints.

Metasurface with circular holes in a square lattice

An often-used design for a resonant metasurface that exhibits a high-Q resonance is a square array of circular holes with a diameter D and a center-to-center pitch Λ^39,76. When certain symmetry conditions are fulfilled, such arrays may support so-called symmetry-protected and/or accidental BICs, depending on the angle of incidence and the exact arrangement of holes. The spectral range at which such a square array supports a BIC is determined by the material and geometrical parameters and it is critical to ensure they are chosen to facilitate both sharp resonances as well as efficient baseline SHG signal. The 3R-MoS₂ flake thickness analysis indicates that there are several optimal thicknesses, which are candidates for metasurface nanofabrication. One may, in principle, aim at the second optimum value, at  ~100–120 nm flake thickness with the strongest SHG signal or the thinner flake with approximately 50% smaller SHG. In the former case, it may be possible that SHG emission will not be able to outcouple efficiently from a thick metasurface. Furthermore, thick samples may present difficulties in achieving precise nanofabrication. This would deteriorate the quality factor of the desired resonances, which are very sensitive to the roughness of the sample⁷⁵. On the other hand, localizing a BIC mode or other high-Q resonances in structures with sub-20 nm thicknesses may be challenging due to an insufficient amount of the material (see Supplementary Fig. S17 for thickness dependence). Balancing the above benefits and limitations, we aim to engineer metasurfaces supporting high-Q optical resonances in multilayer 3R-MoS₂ films of optimal thicknesses – approximately h ~ 20−25 nm, Fig. 2⁷⁷.

Following experimental measurements in Fig. 2c, we aim to generate modes in the 1.3–1.5 eV range, where a thin, unpatterned  ~20 nm flake provides an optimized SHG baseline. With the 3R-MoS₂ thickness set, the remaining parameters, the lattice and hole size, are used to tune the metasurface modes to the desired energy range. Supplementary Fig. S18 demonstrates how the response is red-shifted with increasing lattice and allows us to move the lowest energy modes from approximately 1.3 eV for Λ = 600 nm to 1.6 eV for Λ = 400 nm for thickness 20 nm and hole filling fraction of 0.3. The last parameter is the hole, whose size can be used to predominantly increase the visibility of the mode as well as, to a much lesser degree, provide fine-tuning of the frequency response by  ~0.1 eV, Supplementary Fig. S19. Lastly, the thickness of 3R-MoS₂ will also determine the mode, as its increase will result in a significant red shift of the modes, Supplementary Fig. S17. Hence, precise control of the lattice and thickness is key to setting the desired spectral range, while the hole size is less important in this aspect. However, its value is important in ensuring the proper visibility of the modes.

Symmetry-protected BICs have been realized in square lattices of circular holes in a SiN metasurface at normal incidence⁷⁸. These resonances feature ultra-high quality factors but are challenging to excite from the far field. To effectively excite a BIC, one typically introduces a symmetry break, which can be achieved by illuminating the metasurface at an oblique incidence⁷⁸ or by rotating the meta-atoms within the metasurface⁴⁴. To implement the former, we employ back-focal plane (BFP) imaging and spectroscopy using an oil immersion objective to access the reflection of the metasurface at various incidence angles and polarizations. Upon reaching the normal angle of incidence, the signal related to the BIC mode should vanish from the reflection spectrum, indicating the divergence of the quality factor, as expected for symmetry-protected BICs.

The circular hole metasurface was fabricated using electron beam lithography (EBL) and reactive ion etching (RIE) on a multilayer 3R-MoS₂ flake mechanically exfoliated from a macroscopic crystal (HQ Graphene, additional details are presented in Methods). Figure 3 shows an SEM image of a circular hole array and BFP reflection spectra from this array fabricated in a 27 ± 1-nm thick 3R-MoS₂ flake. The thickness of the flake was measured using AFM.

[See PDF for image]

Fig. 3

Circular hole metasurface.

a Scanning electron microscopy (SEM) image of the circular hole 3R-MoS₂ metasurface with h = 27 ± 1 nm. The scale bar is 1 μm. The inset shows a zoomed-in image of the metasurface. The scale bar is 400 nm. b High-resolution transmission electron microscopy image of one of the circular holes in a 3R-MoS₂ metasurface with a beam direction along the [0001] zone axis. The scale bar is 10 nm. The inset shows the fast Fourier transform (FFT) of the image that exhibits distinct spatial frequency spots corresponding to the crystalline structure of the 3R-MoS₂ and diffuse rings corresponding primarily to the amorphous SiN membrane support. The scale bar is 5 nm⁻¹. The red and green circles identify the FFT spots that correspond to the 3R-MoS₂ {1100} and {1120} families of crystallographic planes, respectively. c Experimental back focal plane reflection spectroscopy of the metasurface for transverse electric (TE) and transverse magnetic (TM) modes with 2 close-ups. The `SHG lens' indicates the effective numerical aperture of the lens during second harmonic generation experiments in air. CLR circular lattice resonance, BIC bound-state-in-the-continuum. Extracted parameters are D_hole = 278 ± 8 nm, Λ_x = 493 ± 7 nm, and Λ_y = 488 ± 8 nm. d Numerically calculated reflection spectra at different angles of reflection for TE and TM modes with similar close-ups. The parameters of the metasurface used in the calculation are: pitch Λ = 500 nm, D = 262.5 nm, h = 23 nm.

The metasurface, presented in the SEM image Fig. 3a, exhibits some shape imperfections and photoresist residues as a result of the non-ideal nanofabrication process and the usage of e-spacer (see Methods and Supplementary Fig. S2). These features may negatively affect the overall performance of the metasurface⁷⁵. The comprehensive statistical analysis of the individual circular hole size (D_hole), shape (ellipticity) and pitch size (Λ_x and Λ_y) of the lattice is described in the Supplementary Fig. S3. Extracted parameters are D_hole = 278 ± 8 nm, Λ_x = 493 ± 7 nm, and Λ_y = 488 ± 8 nm, while the shape of the holes was slightly elliptical. The reasonable resolution, hence error estimation, of the size and shape of the circular hole metasurface from SEM imaging is, conservatively, of the order of 10 nm, which is a substantial value to diminish a high-Q state.

To gain more clarity about the morphology and the atomic structure of the etched holes, we performed high-resolution transmission electron microscopy (HRTEM) and scanning TEM (HRSTEM). Figure 3b shows an HRTEM image of the edge of one of the circular holes on a  ~20-nm thick SiN membrane support. A fast Fourier transformation (FFT) of Fig. 3b is shown in the inset image. The indexed FFT spots reveal the single crystal nature of the 3R-MoS₂ while the diffuse rings primarily correspond to the amorphous SiN membrane support (see Supplementary Fig. S4). The HRTEM image confirms that the 3R-MoS₂ crystallinity extends almost to the edge and there are irregularities in the circular hole edges with a roughness of around 5 nm. Furthermore, there is a  ~10-nm thick region near the hole edge (dark contrast in Fig. 3b) where there is a modification to the morphology and structure. The 3R-MoS₂ crystalline structure continues through this contrast layer, indicating that the contrast is probably residual material on the flake surface near the edge that is left over from the etching process. Additional HRSTEM images (Supplementary Fig. S5) further support the conclusions about the 3R-MoS₂ circular holes.

Figure 3c presents the BFP reflection spectroscopy map of the circular holes metasurface (see Methods and Supplementary Fig. S6 for additional details). The left and right sections correspond to transverse electric (TE) and transverse magnetic (TM) modes, respectively. We observe several narrow dispersive reflection minima corresponding to circular metasurface lattice resonances (CLR) and BIC modes. Specifically, due to the symmetry of the system, for circular nanohole metasurface, the TE and TM polarization modes are expected to align when reaching the normal angle of incidence. The region of interest in terms of angle of incidence is limited to $\sin \theta \approx 0.66$, which corresponds to the total internal reflection angle in the glass-air interface, calculated as $\theta =\arcsin \left(n_{\mathrm{air}}/n_{\mathrm{glass}}\right)\approx 41.5^{\circ }$. The orange and red boxes indicate close-up images of the desired region measured in BFP reflection.

In the case of the TE mode (as shown in the orange box inset in Fig. 3c), two symmetric reflection dips gradually thin down and eventually vanish near normal incidence at photon energies of approximately ℏω ~1.4 eV and 1.45 eV, indicating the appearance of symmetry-protected BICs. In contrast, for the TM mode (indicated in the red box in Fig. 3c), a single flat band exists within the entire range of angles in the inset, indicating that this band is not of a BIC character.

The BFP measurements are complemented by numerical calculations of reflection using the rigorous coupled-wave analysis (RCWA, see Methods for more details), see Fig. 3d. This allows us to extract the optical parameters of the metasurface whose spectrum would match the experimental data. The parameters used in the calculations are Λ = 500 nm, D = 262.5 nm, h = 23 nm, and reasonably close to those extracted for the fabricated sample (see Methods). In addition, the calculations confirm our experimental observations, demonstrating a good agreement in terms of position and dispersion for both TM and TE modes, as is seen from the comparison of Fig. 3c, d. Furthermore, our calculations confirm the CLR and BIC nature of the modes observed at normal incidence at  ~1.4−1.5 eV. To conclude this section, while the reflection spectroscopy of the circular hole 3R-MoS₂ metasurfaces shows promising performance, there remain several aspects to improve, particularly concerning the quality of individual meta-atoms and the precision of their edges.

Quasi-BIC metasurface with triangular holes

As a next step, we fabricate a 3R-MoS₂ metasurface consisting of a square array of triangular holes using anisotropic wet etching of an initially pre-fabricated array of circular holes. This process follows the strategy reported previously for 2H-WS₂ and 2H-MoS₂ (for comparison between the etching of 2H- and 3R-MoS₂, see Supplementary Fig. S7; and for faulty stacked 3R-MoS₂, see Supplementary Fig. S8)^74,79,80 and aims to accomplish two key objectives: (i) achieving ultrasharp edges of the triangular meta-atoms, and (ii) breaking the in-plane C₂ and C₄ rotational symmetry. The former aims at minimizing the scattering losses at the meta-atoms edges, while the latter is necessary to realize the q-BIC concept⁴⁴. To ensure that the metasurface composed of triangular holes offers modes with similar dispersion in the same spectral range as those for the circular hole array, we keep the fill factor (ff) similar to that of the metasurface shown in Fig. 3. Since the material is consumed during wet etching, it is necessary to use a smaller diameter of the initial circular hole to achieve the required triangular hole size. The typical dependence between the side length a_tri of the equilateral triangle hole and the diameter of the initial circular hole D for a given thickness of the 3R-MoS₂ flake (h = 24.5 ± 1.5 nm) shows a linear dependence: a_tri = 2.14D + 72 nm, where D is the hole diameter in nanometers (Supplementary Fig. S9).

It is important to note that material removal occurs in a specific manner: the edges of the triangular hole are oriented along the zigzag (ZZ) edge of the 3R-MoS₂ crystal⁷⁹. Therefore, to achieve a triangular hole array in a perfect square lattice with one side of the triangle parallel to one of the axes of the unit cell, one must know the exact orientation of the ZZ and armchair (AC) directions of the crystal. Moreover, the strict correspondence between lattice orientation and crystalline coordinates plays a pivotal role in SHG experiments, as the χ⁽²⁾ tensor aligns with the crystalline axis, peaking in the AC direction, while optical resonances follow the metasurface lattice orientation. In this study, we aim to optimize both material properties and optical resonances; therefore, it is greatly beneficial that the orientations of nonlinearity and metasurface lattice are aligned.

Figure 4a shows a high-resolution SEM image of a wet-etched metasurface with triangular holes. After the wet etching, the 3R-MoS₂ surface appears to be significantly cleaner than in the case of circular holes, while the edges of the triangular holes are exceptionally sharp (approaching atomic sharpness, due to the self-limiting nature of the wet etching process^74,79). The average size of the triangular hole side lengths (a_tri), which was extracted from the SEM image in Fig. 4a, is a_tri = 373 ± 7.6 nm (Supplementary Fig. S9). Finally, we note that the anisotropic wet etching is able to reveal stacking faults that occasionally occur in 3R-MoS₂ multilayers (Supplementary Fig. S8); such faulty metasurfaces were excluded from reflection and SHG experiments. Furthermore, the apparent absence of defects in the triangular metasurface shown in Fig. 4 suggests that it is entirely free from stacking faults and 2H inclusions. Since the likelihood of stacking faults increases with sample thickness, this observation further supports the use of ultrathin metasurfaces with h ~ 20–25 nm.

[See PDF for image]

Fig. 4

Triangular hole metasurface.

a High-resolution scanning electron microscopy (SEM) image of the h = 24.5 ± 1.5 nm thick metasurface. The scale bar is 1 μm. Inset is a close-up SEM image of the metasurface. The scale bar is 400 nm. b(i) High-resolution transmission electron microscopy image of a vertex of one of the triangular holes in a 3R-MoS₂ metasurface with a beam direction along the [0001] zone axis. The scale bar is 5 nm. The inset shows the fast Fourier transform (FFT) of the image that exhibits distinct spatial frequency spots corresponding to the crystalline structure of the 3R-MoS₂. The scale bar is 5 nm⁻¹. The red and green circles identify the FFT spots that correspond to the 3R-MoS₂ {1100} and {1120} families of crystallographic planes, respectively. The orientation between the {1100} family of spots in the FFT and the edges of the triangle in the image confirm the zigzag nature of the etched edges. (ii) High-resolution scanning transmission electron microscopy image of a representative triangular hole edge. The overlays show possible zigzag and armchair directions, confirming the edge is aligned with the zigzag direction. The scale bar is 2 nm. c Experimental back focal plane reflection spectroscopy of the metasurface for transverse electric (TE) and transverse magnetic (TM) modes with 2 close-ups. TLR triangular lattice resonance, qBIC quasi-bound-state-in-the-continuum. d Numerically calculated reflection spectra at different angles of reflection for TE and TM modes with similar close-ups using following parameters: a_tri = 375 nm, pitch size Λ = 500 nm, h = 20 nm.

HRTEM and HRSTEM images of triangular holes, as shown in shown in Fig. 4b and Supplementary Figs. S10–S14 confirm that the triangular hole edges have better crystallinity and are sharper compared to the circular hole edges. The HRTEM image of a triangle vertex in Fig. 4b(i) and the HRSTEM image of a triangle edge in Fig. 4b(ii) show that the triangle vertex is exceptionally sharp and that the edges are nearly atomically sharp with a roughness ≤1–2 nm. The orientation between the edges in the HRTEM image in Fig. 4b(i) and the indexed FFT spots in the inset reveal that the triangle edges are aligned with the zigzag direction. This is also confirmed by the 3R-MoS₂ atomic structure that is resolved in the HRSTEM image in Fig. 4b(ii), as indicated by the overlay markers. Additional HRTEM and HRSTEM images (Supplementary Figs. S10–S14) further support the conclusions about the 3R-MoS₂ triangular holes.

As with the circular hole array, we measure reflection spectra in BFP for the triangular hole metasurface. Unlike the former case, the unit cell with an equilateral triangular hole lacks in-plane C₂ and C₄ symmetries, therefore a metasurface based on such a cell is a promising candidate to allow the existence of the q-BIC state at normal incidence from the far-field excitation⁴⁴. The collected experimental BFP spectra are plotted in Fig. 4c for TE and TM polarizations with the same orientation as for the circular hole metasurface. Overall, the mode dispersion for the triangular metasurface shows similar behavior to the circular one, however, with a few exceptions. Specifically, unlike the circular hole metasurface, in this case, a closer look reveals the presence of a reflection dip for the lower order mode even at normal incidence at 1.47 eV, indicating the presence of q-BIC. This normal incidence q-BIC is confirmed in the calculated reflection spectra for a metasurface with the following parameters: a_tri = 375 nm, lattice Λ = 500 nm, h = 20 nm (Fig. 4d), which closely match the parameters extracted experimentally: a_tri = 373 ± 7.6 nm, lattice pitches Λ_x = 501 ± 7.3 nm, Λ_y = 507 ± 7.4 nm, and flake thickness h = 24.5 ± 1.5 nm (see Methods). Both experimental and theoretical results prove that the wet-etching fabrication is a powerful tool to break the in-plane symmetry of the meta-atom to access the q-BIC modes alongside reaching exceptional sharpness of the meta-atoms’ edges.

SHG enhancement in circular and triangular hole metasurfaces

After careful analysis of SHG in bare 3R-MoS₂ flakes and BFP reflection of metasurfaces, which identified the existence of promising lattice resonance and optical BIC (q-BIC) modes in the 1.4–1.6 eV range, we proceeded to study the SHG in those metasurfaces. Figure 5a shows SHG enhancement spectra of the circular hole metasurface for a linearly polarized excitation along the y-axis (specified in the inset). Circular points are experimental SHG enhancement data, namely the ratio of the SHG intensity of the metasurface over the SHG intensity of the host flake at a given pump wavelength. We introduce this figure of merit – SHG enhancement – to quantify how the high-Q modes can amplify the SHG signal of the metasurface compared to the SHG signal of the optimized host flakes. However, since χ⁽²⁾ of 3R-MoS₂ depends on the pump wavelength, direct comparison between metasurfaces with resonances occurring at different wavelengths is not straightforward. One can see that the circular hole metasurface features a clear resonant enhancement around 850 nm, which corresponds to a pump energy of 1.46 eV – the same energy as the supported TM mode observed in the BFP spectra. The experimental SHG enhancement reaches  ~200-fold, while the numerically calculated one is almost an order of magnitude larger at the same pump wavelength (gray shaded area in Fig. 5a). The numerically calculated SHG of the metasurface was carried out using the same parameters as the calculated reflection spectra in Fig. 3c. The SHG of the host flake was calculated using the same thickness as the metasurface (for additional details about SHG calculations, see Methods). The mismatch in the magnitude of SHG enhancement between the calculation and experiment can be attributed to the surface and shape roughness of the circular holes, as one can see from the AFM image shown in the inset of Fig. 5a, and from the SEM image in Fig. 3a.

[See PDF for image]

Fig. 5

Second harmonic generation (SHG) of circular and triangular hole metasurfaces.

a Circular hole metasurface SHG enhancement, SHG_meta/SHG_flake, vs. pump wavelength λ_pump. The pump polarization is linear along the y-axis. Detection has no analyzer. Circles correspond to experimental data, while the shaded area to numerical simulations. The inset depicts an atomic force microscopy (AFM) scan of the metasurface. b Triangular hole metasurface SHG enhancement, SHG_meta/SHG_flake, vs. the pump wavelength. The pump polarization is linear along y-axis. Detection has no analyzer. Triangles correspond to experimental data, while the shaded area to numerical calculations in Comsol. The inset depicts an AFM scan of the metasurface. The error bars in (a) and (b) are less than the marker sizes. c Numerical (calc) - (i, iii) - and experimental (exp) - (ii,iv) polar plots of the SHG intensity (i, ii) and SHG enhancement (ratio of metasurface-to-flake SHG; iii,iv)) for the circular hole metasurface. SHG intensity (i,ii) is plotted in a linear scale. SHG enhancement (ratio; iii, iv) is plotted in a log-scale. The shaded background is the AFM image of the metasurface aligned with the polarization angles. d Numerical (calc) - (i, iii) and experimental (exp) - (ii, iv) polar plots of the SHG intensity (i, ii) and SHG enhancement (ratio of metasurface-to-flake SHG; iii, iv) for the triangular hole metasurface. SHG intensity (i, ii) is plotted in a linear scale. SHG enhancement (ratio; iii, iv) is plotted in a log-scale. The shaded background is the AFM image of the metasurface aligned with the polarization angles. Paramaters of the numerical calculation are the same as in Figs. 3 and 4.

Similar experiments and numerical analysis were carried out for the triangular hole metasurface. Figure 5b shows SHG enhancement of the triangular hole array with the incident polarization along the y-axis (along the height of the triangular hole and parallel to one of the axes of the unit cell). Similarly to the circular hole case, the spectral match is nearly perfect between the experimental SHG enhancement and the numerical one with the parameters extracted from Fig. 4d, peaking at roughly λ_pump = 810 nm, which matches with the mode appearing in BFP around 1.53 eV in Fig. 4c. However, in the case of the triangular hole metasurface, numerical and experimental data match not only spectrally, but also in the magnitude of the SHG enhancement –  ~800-fold experimentally and  ~1400-fold for the numerical study. It is worth mentioning, that the numerically perfect circular and triangular hole arrays perform similarly in terms of the peak SHG enhancement (~1500-fold), however, experimentally we obtain substantially larger values for the triangular case –  ~800-fold vs.  ~200-fold. We associate these results with the sharpness of the edges in triangular nanoholes and the cleanliness of the sample surface, as one can see in the AFM scan of the triangular hole sample with ultrasharp edges and a residual-free surface (Fig. 5d inset).

After optimizing the pump wavelength for the circular hole metasurface at 850 nm, we perform polarization-resolved SHG measurements at this excitation wavelength (see Methods). Figure 5c shows experimental data of the SHG intensity of the circular hole metasurface (blue area) as a function of the angle ϕ characterizing the rotation of pump polarization with respect to the x-axis, and SHG of the nearby bare host flake (red line) at 850 nm. The black circles are the ratio between the SHG of the metasurface and the host flake indicating how efficiently the metasurface can enhance the SHG intensity at a given angle. In Fig. 5c one observes at ϕ = 90°, which corresponds to the maximal SHG in both the host flake and the lattice resonance of the metasurface (along AC axis of the flake), the enhancement factor, although smaller than maximum values, exceeds 3 orders of magnitude.

The Fig. 5c(i) shows the calculated SHG intensity of the metasurface and the host flake at λ_pump = 850 nm overlaid with a schematic of the metasurface to illustrate the alignment of polarization-resolved SHG signal with respect to the metasurface axes. The data indicate perfect six-fold symmetry of the flake (red line) and 3 pairs of 2-fold symmetric features in the SHG of the metasurface with one dominant direction (blue line, dominant along ϕ = 90° direction). Since the total SHG intensity is a product of the field enhancement and the nonlinearity of the material system, it is evident that the vertical direction is dominant for SHG in terms of both the AC direction of the flake and the lattice resonance of the metasurface. However, the appearance of non-zero SHG along the directions at which the host flake has theoretically vanishing SHG (ZZ-axis) provides substantial SHG enhancement (up to 8 orders of magnitude in calculations), which is plotted in Fig. 5c(iii) (black line).

The Fig. 5c(ii) shows experimental polarization-resolved SHG for circular hole metasurface in polar coordinates. The experimental data exhibits substantially less pronounced subpeaks (blue line, at ~30° and 330°), however, they appear along approximately the same direction as in the numerical calculations. Such a discrepancy could be caused by (i) a slight misalignment of the crystalline axes and the metasurface unit cell and/or (ii) the metasurface not supporting high-Q modes because of the shape and surface imperfections. The latter hypothesis is also confirmed by the modest SHG enhancement shown in the Fig. 5c(iv) (black line), where the maximum enhancement does not reach 3 orders of magnitude along the ZZ direction where the reference flake exhibits minimum SHG.

The polarization-resolved analysis of the SHG enhancement in triangular metasurface is shown in Fig. 5d. The panel (i) shows a numerical calculation of the SHG intensity of the triangular hole metasurface (blue line) and the host flake (red line) as a function of ϕ. In this case, the SHG intensity of the metasurface exhibits a single 2-fold symmetry feature along the triangular hole array lattice and the AC axis of the host flake – y-axis. This is in agreement with the experimental result, shown in panel (ii), however, the experimental polar plot exhibits a slight tilt in the peak position (by several degrees in a clockwise direction). This tilt is also observed in AFM (inset Fig. 5b) and SEM (Fig. 4a) images of the triangles and is likely caused by a slight misalignment of the crystalline axes of the 3R-MoS₂ host flake with respect to the metasurface’s unit cell, which in turn induces a slight symmetry break. Finally, the calculated SHG enhancement in the Fig. 5d(iii) (black line) is in good agreement with the experimental data in panel (iv) (black line) supporting the dominant direction of the enhancement along the ZZ direction of the flake and the experimentally measured SHG enhancement factor along ϕ ≈ 60° direction exceeds 10³. In summary, the triangular nanohole metasurface presented here demonstrates performance that is considerably improved in comparison to its circular counterpart, making it a promising candidate for future nonlinear nanophotonics applications.

To additionally demonstrate the final point and benchmark our metasurfaces against established standards, we directly compared both circular and triangular metasurfaces to SHG from a monolayer MoS₂. Our results demonstrate that the SHG signal is enhanced ~33,000- and 133,000-fold for the circular and triangular metasurfaces, respectively, relative to monolayer MoS₂, values which are achieved despite a moderate increase in the overall sample thickness. Furthermore, the absolute conversion efficiencies, defined as η_∘,△ = P_SHG/P_pump, are on the order of η_∘ = 0.02% and η_△ = 0.08% for circular and triangular metasurfaces, respectively. Note that these conversion efficiencies are estimated at a moderate P_pump = 0.3 mW, at 800–850 nm pump wavelength range, repetition rate of 80 MHz, and pulse duration of  ~100 fs (see Methods for further details).

Conclusions

In this work, we used the large refractive index and non-centrosymmetric features of 3R-MoS₂ to fabricate efficient ultrathin (h = 20–25 nm) nonlinear metasurfaces. First, we analyzed the thickness and wavelength dispersion of pristine 3R-MoS₂ flakes for SHG in the h = 0–400 nm thicknesses and λ_pump = 800–1040 nm pump wavelength (1.2–1.55 eV) range. Second, we demonstrated the wet-etching fabrication of atomically sharp triangular hole metasurfaces with broken in-plane C₂ and C₄ symmetries, supporting q-BIC resonances. The atomic sharpness of the meta-atoms’ edges alongside the stacking-fault-free nature of our metasurfaces leads to a substantial enhancement of the SHG signal, exceeding three orders of magnitude at specific resonant wavelengths and linear polarization conditions. This enhancement is achieved with respect to the substantial SHG signal of the optimally thick host flake (corresponding to the first optimum thickness range, h = 20−25 nm, Fig. 2).

We note that the metasurface resonances fall within the 800–900 nm range, overlapping with the fundamental wavelength range of our tunable femtosecond laser. Consequently, in this study, we focus on optical resonances that overlap with the fundamental wavelength rather than the second harmonic. Nevertheless, our platform has the potential to support more complex configurations, where optical resonances could be, in principle, tuned to overlap with either the fundamental or second-harmonic frequency, or even both simultaneously.

Finally, our fabrication technique offers precise control over the photonic lattice, which is rigidly linked with the crystalline axes of 3R-MoS₂. The relative orientation between the crystal and metasurface axes introduces an additional degree of freedom for symmetry control, potentially enabling chiral harmonic generation, as predicted recently⁸¹. Furthermore, 3R-MoS₂ metasurfaces may enable active tuning of the nonlinear emission response via electro-optical modulation, carrier injection, and piezoelectric effects^66,67. Thus, our work opens the door to nonlinear nanophotonics in high-quality van der Waals metasurfaces, at the interplay of atomically sharp edges, high refractive index, and extreme nonlinearities.

Methods

Sample fabrication

The nanofabrication steps were carried out in the Myfab Nanofabrication Laboratory, MC2 Chalmers. A high-quality 3R-MoS₂ crystal was purchased from HQ-Graphene and mechanically exfoliated into multilayer flakes. Multilayer flakes with the desired thickness were selected and transferred to 1-inch  × 1-inch glass substrates (0.17 mm thickness) with the help of PDMS stamps (Gel-Pak, USA). PMMA was used as a positive resist for electron beam lithography (EBL). A Raith EBPG 5200 (Germany) system operating at 100 kV accelerating voltage was used for the direct-writing EBL. The PMMA mask was nanopatterned with holes using 10 nA beam current. Oxford Plasmalab 100 system (RIE/ICP) was used with CHF₃/Ar to perform the dry etching step after EBL to transfer the pattern through the mask to the 3R-MoS₂ flakes. The leftovers of the resist were removed by a remover and, finally, the samples were washed with deionized water and gently blow-dried. For the circular hole arrays, the nanofabrication steps ended here. However, for the triangular hole arrays, an additional step was required to transform the initial circular holes into triangular shapes. This conversion was achieved through anisotropic wet etching using an aqueous solution of hydrogen peroxide and ammonia, as detailed in ref. ⁷⁴. The nanofabricated devices were then stored in a cleanroom environment until further experiments.

For the preparation of TEM samples, circular and triangular hole arrays were first prepared on SiO₂/Si surfaces using the same procedure as described above. The patterned flakes were picked up using PDMS stamps and subsequently transferred onto silicon nitride (SiN) membrane TEM grids (simpore.com and www.tedpella.com) using an all-dry transfer process to ensure high quality. Immediately before the transfer of flakes, the SiN membranes were plasma cleaned in a Fischione plasma cleaner in a mixture of 25% oxygen and 75% argon for 30s to remove organic contaminants from the surface of the SiN membrane. This improved the adhesion between the membrane and the flakes during flake transfer. The circular patterned flakes were transferred onto continuous 20-nm-thick SiN membranes. The triangular patterned flakes were transferred onto 200-nm-thick holey SiN membranes featuring pre-patterned holes with a diameter of 600 nm. Some of the triangular etched patterns were positioned directly over these holes to enable direct observation of the atomic structure of the triangular etched edges of 3R-MoS₂ without interference from the underlying SiN substrate. The samples were promptly loaded into the microscope immediately after the TEM sample preparation.

Atomic force microscopy

Atomic Force Microscopy (AFM) imaging was performed using a Bruker ICON AFM system equipped with a Nanoscope 5 controller, operating in tapping mode. Topography measurements were carried out using Bruker RTESP-300 AFM probes.

Scanning electron microscopy

SEM imaging of the circular and triangular hole metasurfaces was performed at Chalmers Materials Analysis Laboratory using an Ultra 55 microscope (Carl Zeiss). E-spacer 300Z was spin-coated onto the samples as a conductive layer to reduce charging effects, as the samples were fabricated on a non-conductive glass substrate. An acceleration voltage of 1.7 kV was used for imaging both samples.

The SEM images were analyzed using ImageJ software, where a monochromatic threshold filter was employed to distinguish the circular/triangular holes from the basal plane of the flake. The surface area S of each circular/triangular hole was analyzed by ImageJ (450 holes in total). For the case of circular holes, the cross sections were fitted by an ellipsoid to obtain the major and minor axes, as well as the tilt angle between the major axis and the x-axis. The diameter distribution, which is plotted in Supplementary Fig. S3, is derived from the surface area of each hole (D = 2$\sqrt{S/\pi }$) assuming perfect circles. The analysis reveals pitch sizes of Λ_x = 493 ± 7.4 nm, Λ_y = 488 ± 8.3 nm, the mean diameter of the circular holes is ⟨D⟩ = 278 ± 8 nm, and the fill factor ff = 0.25 ± 0.02, which is defined as ff = S/Λ², i.e. the ratio of the area of a circular hole to the area of the square unit cell.

A similar analysis for triangular metasurface reveals the surface area S_△ of each of the 450 holes, which, assuming equilateral triangles, results in the side lengths as $a_{\mathrm{tri}}=\sqrt{2S_{△}/\sin \left(\pi /3\right)}$. The mean side length is ⟨a_tri⟩ = 373 ± 7.6 nm. It is worth mentioning that the relative standard deviation of the triangular holes’ size is smaller than that of circular holes, being 2% and 3%, respectively. This is illustrated in Supplementary Fig. S9.

We note that the error bars appearing in SEM image processing are dominated by systematic errors caused by the image resolution limitations of low-magnification images (~14 nm/pixel). Higher resolution SEM images, which contain enough nanoholes to perform meaningful statistics, were challenging to record in this case, due to significant charging effects, occurring because the samples resided on non-conductive SiO₂ substrates. E-spacer was used to reduce the charging effects, however, it was not possible to fully avoid the issue. For the triangular nanoholes, several higher-resolution SEM and TEM images are shown in the SI (Supplementary Figs. S7, S10–S14), demonstrating that the edges of these holes are exceptionally sharp.

Transmission electron microscopy

We performed HRTEM and HRSTEM imaging to analyze the morphology and the atomic structure of the circular and the triangular etched 3R-MoS₂ flakes. HRTEM experiments were conducted on a Thermo Fisher Titan microscope at an accelerating voltage of 300 kV. HRSTEM experiments were conducted on a JEOL Mono NEO ARM 200F at an accelerating voltage of 200 kV. The instrument is equipped with a Schottky field emission gun, double Wien filter monochromator, probe aberration corrector, image aberration corrector, and Gatan Imaging Filter continuum HR spectrometer. HRSTEM images were restored to reduce the effects of image noise and distortions using in-house developed python scripts that utilize deep convolution neural networks⁸². Supplementary Fig. S12 shows a visual comparison of a raw and restored image. Double tilt holders were used to align the incident electron beam with the MoS₂ [0001] crystallographic zone axis.

Optical reflection spectroscopy measurements

Linear optical reflection spectroscopy was carried out using a back-focal plane reflection optical setup. The setup consists of an inverted microscope (Nikon Eclipse TE2000-E) with an oil immersion objective NA = 1.49 (Nikon CFI Apo TIRF oil, MRD01691, 60×). A laser-driven broadband light source (LDLS, EQ-99FC) was used to illuminate the samples from the glass substrate’s side with immersion oil (Nikon NF, n_oil = 1.515). The immersion oil allows access to high NA with a theoretical limit of ${\theta }_{\max }=\arcsin \left(\mathrm{NA}/n_{\mathrm{oil}}\right)$. The Fourier plane of the objective was imaged using a Bertrand lens and its images were recorded using a digital color camera (Nikon D300S). The BFP reflection spectra were recorded by a spectrometer (Andor Shamrock SR-500i, equipped with a CCD detector Andor Newton 920), while the fiber bundle consisting of 19 individual fibers allows us to collect 19 individual reflection spectra simultaneously (the schematic of the setup is shown in the SI, Supplementary Fig. S6).

Second-harmonic generation spectroscopy measurements

SHG experiments performed on bare 3R-MoS₂ flakes, presented in Fig. 2, were carried out using an air NA = 0.95 objective with a correction ring (Nikon CFI Plan Apochromat Lambda D, 40×) to compensate for distortion occurring by passing through the glass layer. The SHG of the metasurfaces was measured using an air NA = 0.1 objective (Nikon CFI Plan Achro, 4×). The low NA objective was chosen to mimic a near-normal incidence excitation, where angular dispersion of the metasurface resonances is expected to be weak (see BFP reflection data).

Polarization-resolved SHG data presented in Fig. 5a, b was obtained using a linearly-polarized laser beam and without an analyzer installed in the signal collection arm. Polarization-resolved SHG data presented in Fig. 5c, d was obtained using a linear polarization rotation with a λ/2 plate with 5° steps. The SHG beam of the reflected beam passed the same sequence of λ/2 and linear polarizer as the incident light before being detected by the APD, therefore leaving only the component parallel to the incident linear polarization.

In all SHG experiments, a tunable, 690–1040 nm, Ti : sapphire femtosecond laser (MaiTai HP-Newport Spectra-Physics) with a  ~100 fs pulse duration and 80 MHz repetition rate was used. The SHG signal was collected through a fiber and analyzed by a spectrometer (Andor Shamrock SR-500i, equipped with a CCD detector Andor Newton 920) for single spectra, and with an avalanche photodiode (APD, IDQ id100) for mapping using piezo-stage (Mad City Labs, MCL Nano-LP200) for precise positioning of the sample. The explicit SHG spectra, alongside SHG intensity vs. input laser power demonstrating the slope = 2 in the log-log scale (and thus proving the second-harmonic nature of the recorded signal) are presented in the SI (Supplementary Fig. S1).

Benchmarking metasurfaces vs. monolayer MoS₂ and absolute SHG conversion efficiencies

Our results reveal that the SHG signal is enhanced by approximately 33,000- and 133,000-fold for the circular and triangular metasurfaces, respectively, relative to monolayer MoS₂. Part of this enhancement arises from the thickness effect (propagation length, h, inside the nonlinear medium), while the rest is due to the optical resonances in the metasurfaces. Specifically, the SHG of the metasurfaces was measured using an air 4×, NA = 0.1 objective in reflection configuration (thus, both the pumping and SHG collection were performed through the same objective), and for benchmarking, we focus on the 800–900 nm range of the pumping laser with a  ~100 fs pulse duration and 80 MHz repetition rate. The conversion efficiencies, η, for circular and triangular metasurfaces were on the order of 0.02% and 0.08%, respectively, at moderate P_pump ~ 0.3 mW (corresponding to  ~5 × 10⁷ W/cm² peak power density). P_SHG and P_pump were measured by an APD using the SHG signal and attenuated laser light, respectively, with the detector’s quantum efficiency at the second-harmonic and fundamental wavelengths taken into account (Supplementary Fig. S20, APD calibration). The damage threshold of the 3R-MoS₂ flakes are on the order of 50–70 mW for similar excitation geometry⁶⁸. Since the conversion efficiency η ∝ P_pump, η may be realistically enhanced above  ~0.08% in triangular metasurfaces before reaching the damage threshold. Such high conversion efficiencies in ultrathin metasurfaces (20–25-nm thickness) indicate an exceptionally efficient nonlinear process, approaching that of more traditional bulk nonlinear materials.

Calculation of linear response of metasurfaces

Reflection spectra of metasurfaces were calculated using the rigorous coupled wave analysis (RCWA) assuming a square periodic lattice. The geometrical parameters were taken from fabrication designs, structural characterization, and SEM images and then varied in a narrow parameter range to obtain good matching to experimental spectra. The permittivity of the semi-infinite glass substrate was assumed to be 2.25 and that of MoS₂ was taken from ellipsometric measurements⁵⁴. The resolution of the patterned MoS₂ layer was 2.5 nm and we took up to 250 terms for the plane wave expansion of the fields, which ensured converged results. For circular holes, we assumed azimuthal incidence along one of the principal axes of the square lattice. For the triangular pattern, we considered incidence both along the base as well as the height of the triangle.

Linear response calculations reveal resonance quality factors of approximately 100 in the 800–900 nm range (Figs. 3d and 4d, as well as Supplementary Figs. S17–S19). While these values are relatively modest, they are well-matched to the pulse bandwidth of our femtosecond excitation laser used for SHG experiments. Specifically, for our  ~100 fs pulse duration, optimal SHG is achieved at moderate quality factors around 100, which our metasurfaces realistically provide.

Calculation of 3R-MoS₂ bare film SHG

The SHG spectra of 3R-MoS₂ flakes at normal illumination and collection were calculated using the transfer matrix method. The material parameters were assumed as above, while the nonlinear susceptibility was taken from electronic structure calculations reported in Zograf et al.⁶⁸. We calculated the electric field profiles in the MoS₂ layer at the fundamental and second harmonic frequencies. The fundamental mode field profile was used to obtain the macroscopic nonlinear polarization. The nonlinear polarization was then multiplied by the mode profile at the second harmonic and an overlap integral was evaluated to yield the SHG signal at normal emission⁸³.

Calculation of metasurface SHG

Due to D_3h symmetry of the 3R-MoS₂, the total SHG nonlinear polarization can be estimated through the following relation^15,84,85:

2

$\mathbf{P}\left(2\omega \right)=\left[\begin{array}{c}\hfill P_x\hfill \\ \hfill P_y\hfill \\ \hfill P_z\hfill \end{array}\right]\propto \left[\begin{array}{c}\hfill -2\mid {\chi }^{\left(2\right)}\mid E_x\left(\omega \right)E_y\left(\omega \right)\hfill \\ \hfill \mid {\chi }^{\left(2\right)}\mid \left(E_y{\left(\omega \right)}^2-E_x{\left(\omega \right)}^2\right)\hfill \\ \hfill 0\hfill \end{array}\right],$where $\mid {\chi }^{\left(2\right)}\mid ={\chi }_{yyy}^{\left(2\right)}=-{\chi }_{yxx}^{\left(2\right)}=-{\chi }_{xxy}^{\left(2\right)}=-{\chi }_{xyx}^{\left(2\right)}$ – being the dominant tensor components due to symmetry constraints. The χ⁽²⁾ dispersion was taken from ref. ⁶⁸. To reproduce results for polarization-resolved SHG polar plots for normal incidence of the pump along z-axis with just the parallel component of the SHG signal passing through the analyzer, we multiply the resulting SHG signal by the Jones matrix for a given linear polarization angle ϕ:

3

${\mathbf{P}}_{\varphi }\left(2\omega \right)=\left[\begin{array}{cc}\hfill {\cos }^2\left(\varphi \right)\hfill & \hfill \cos \left(\varphi \right)\sin \left(\varphi \right)\hfill \\ \hfill \cos \left(\varphi \right)\sin \left(\varphi \right)\hfill & \hfill {\sin }^2\left(\varphi \right)\hfill \end{array}\right]\left[\begin{array}{c}\hfill P_x\hfill \\ \hfill P_y\hfill \end{array}\right],$where P_ϕ(2ω) indicates the SHG nonlinear polarization along angle ϕ.

The SHG calculation of the periodic arrays (metasurfaces) was carried out in Comsol Multiphysics 6.2. The computational domain consisted of a rectangular box with a square cross-section along xy-plane. The size of the square section is defined by the pitch size of the metasurface Λ. The excitation port was placed on top of the computational domain and the excitation wave was set to propagate along the z-axis. The linear polarization of the excitation field was $E_x=E_0\cos \varphi$ and $E_y=E_0\sin \varphi$, where E₀ is the incident light intensity. The adjacent domain walls perpendicular to the xy plane have a Flouqet periodicity for the k-vector. The computational steps consisted of two physical problems: (i) – fundamental field and (ii) – second harmonic. In step (i), the problem of the electromagnetic field distribution in the metasurface at 800–1040 nm pump wavelengths is solved. The values of the E_x(ω) and E_y(ω) at the fundamental wavelength define the nonlinear polarization used in the second step (see Eq. (2)). Then, in step (ii), the electromagnetic wave distribution problem is calculated at the second-harmonic wavelength. The total SHG intensity collected during the experiment is defined by the Poynting vector penetrating the upper surface of the computational domain. For polarization-resolved SHG calculations and SHG enhancement of the parallel component, we use the same computational model, however, the emitted SHG is then analyzed by taking into account the Jones matrix as described in Eq. (3).

Acknowledgements

G.Z. and B.K. are grateful to Dr. Adriana Canales for assistance with the optical back-focal plane setup. G.Z. thanks Ryo Mizuta Graphics for the optical components library for optical setups visualization. G.Z., B.K., and T.O.S. acknowledge funding from the Swedish Research Council (VR project, grant No. 2022-03347), Chalmers Area of Advance Nano, 2D-TECH VINNOVA competence center (Ref. 2019-00068), and Olle Engkvist foundation (grant No. 211-0063). T.O.S., A.B.Y., A.R., and E.O. acknowledge funding from the Knut and Alice Wallenberg Foundation (KAW, grant No. 2019.0140). M.B. and T.J.A. acknowledge support from the Polish National Science Center via the project 2019/34/E/ST3/00359. A.B.Y. acknowledges funding from the Swedish Research Council (VR starting, grant No. 2020-04986). W.W. acknowledges support from the Knut and Alice Wallenberg Foundation (KAW fellow, grant No. 2019.0231). Calculations were partially done at the ICM, UW (#GC84-51). This work was performed in part at the Chalmers Material Analysis Laboratory (CMAL) and Myfab Chalmers.

Author contributions

G.Z., B.K., W.W., and T.O.S. conceived the idea. G.Z. with the help of B.K. measured all linear and non-linear optical responses of the samples. B.K. with the help of A.P. fabricated all the samples. B.K. performed SEM imaging. A.A. provided AFM imaging of all samples. A.B.Y. and A.R. performed the HRTEM and HRSTEM imaging and analysis. A.B.Y., A.R., and E.O. wrote the TEM-related text. E.O. supervised the TEM-related part of the project. T.J.A. with the help of M.B. and G.Z. investigated numerically linear and nonlinear optical responses of the samples. T.J.A. and T.O.S supervised the theoretical and experimental parts of the project, respectively. G.Z., B.K., and T.O.S wrote the manuscript with the contributions of all co-authors.

Peer review

Peer review information

This manuscript has been previously reviewed at another journal that is not operating a transparent peer review scheme. The manuscript was considered suitable for publication without further review at Communications Physics.

Funding

Open access funding provided by Chalmers University of Technology.

Data availability

The principal data supporting the findings of this study are available within the article and its Supplementary Information files, from the corresponding authors on reasonable request or via Zenodo at https://doi.org/10.5281/zenodo.14194107⁸⁶.

Code availability

Calculations which produce data used in this study were done using Comsol and RCWA. Model files for Comsol version 6.2 and the python script for RCWA calculations are available via Zenodo at https://doi.org/10.5281/zenodo.14194107⁸⁶ or from the corresponding authors upon reasonable request.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s42005-025-02194-y.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Manzeli, S; Ovchinnikov, D; Pasquier, D; Yazyev, OV; Kis, A. 2D transition metal dichalcogenides. Nat. Rev. Mater.; 2017; 2, pp. 1-15. [DOI: https://dx.doi.org/10.1038/natrevmats.2017.33]

2. Mak, KF; Lee, C; Hone, J; Shan, J; Heinz, TF. Atomically thin MoS₂: a new direct-gap semiconductor. Phys. Rev. Lett.; 2010; 105, 136805.2010PhRvL.105m6805M [DOI: https://dx.doi.org/10.1103/PhysRevLett.105.136805]

3. Splendiani, A et al. Emerging photoluminescence in monolayer MoS₂. Nano Lett.; 2010; 10, pp. 1271-1275.2010NanoL.10.1271S [DOI: https://dx.doi.org/10.1021/nl903868w]

4. Wang, G et al. Colloquium: excitons in atomically thin transition metal dichalcogenides. Rev. Mod. Phys.; 2018; 90, 021001.2018RvMP..90b1001W3833231 [DOI: https://dx.doi.org/10.1103/RevModPhys.90.021001]

5. Khan, AR et al. Optical harmonic generation in 2D materials. Adv. Funct. Mater.; 2022; 32, 2105259. [DOI: https://dx.doi.org/10.1002/adfm.202105259]

6. Zhou, X et al. Strong second-harmonic generation in atomic layered GaSe. J. Am. Chem. Soc.; 2015; 137, pp. 7994-7997. [DOI: https://dx.doi.org/10.1021/jacs.5b04305]

7. Yao, K et al. Continuous wave sum frequency generation and imaging of monolayer and heterobilayer two-dimensional semiconductors. ACS Nano; 2019; 14, pp. 708-714. [DOI: https://dx.doi.org/10.1021/acsnano.9b07555]

8. Wang, Y et al. Difference frequency generation in monolayer MoS₂. Nanoscale; 2020; 12, pp. 19638-19643. [DOI: https://dx.doi.org/10.1039/D0NR01994A]

9. Trovatello, C et al. Optical parametric amplification by monolayer transition metal dichalcogenides. Nat. Photon.; 2021; 15, pp. 6-10.2021NaPho.15..6T [DOI: https://dx.doi.org/10.1038/s41566-020-00728-0]

10. Dinparasti Saleh, H et al. Towards spontaneous parametric down conversion from monolayer MoS₂. Sci. Rep.; 2018; 8, 2018NatSR..8.3862D [DOI: https://dx.doi.org/10.1038/s41598-018-22270-4] 3862.

11. Guo, Q et al. Ultrathin quantum light source with van der Waals NbOCl₂ crystal. Nature; 2023; 613, pp. 53-59.2023Natur.613..53G [DOI: https://dx.doi.org/10.1038/s41586-022-05393-7]

12. Sun, Z; Martinez, A; Wang, F. Optical modulators with 2D layered materials. Nat. Photon.; 2016; 10, pp. 227-238.2016NaPho.10.227S [DOI: https://dx.doi.org/10.1038/nphoton.2016.15]

13. Kumar, N et al. Second harmonic microscopy of monolayer MoS₂. Phys. Rev. B; 2013; 87, 161403.2013PhRvB.87p1403K [DOI: https://dx.doi.org/10.1103/PhysRevB.87.161403]

14. Marini, L; Helt, L; Lu, Y; Eggleton, BJ; Palomba, S. Constraints on downconversion in atomically thick films. J. Opt. Soc. Am. B; 2018; 35, pp. 672-679.2018JOSAB.35.672M [DOI: https://dx.doi.org/10.1364/JOSAB.35.000672]

15. Boyd, R. W., Gaeta, A. L. & Giese, E. Nonlinear optics. In Springer handbook of atomic, molecular, and optical physics, 2nd edn., 1097–1110 (Springer, 2023).

16. Ou, Z; Pereira, S; Polzik, E; Kimble, HJ. 85% efficiency for CW frequency doubling from 1.08 to 0.54 μm. Opt. Lett.; 1992; 17, pp. 640-642.1992OptL..17.640O [DOI: https://dx.doi.org/10.1364/OL.17.000640]

17. Bernhardt, N et al. Quasi-BIC resonant enhancement of second-harmonic generation in WS₂ monolayers. Nano Lett.; 2020; 20, pp. 5309-5314.2020NanoL.20.5309B [DOI: https://dx.doi.org/10.1021/acs.nanolett.0c01603]

18. Wang, T; Li, Z; Zhang, X. Improved generation of correlated photon pairs from monolayer WS₂ based on bound states in the continuum. Photon. Res.; 2019; 7, pp. 341-350. [DOI: https://dx.doi.org/10.1364/PRJ.7.000341]

19. Kravtsov, V et al. Nonlinear polaritons in a monolayer semiconductor coupled to optical bound states in the continuum. Light Sci. Appl.; 2020; 9, 2020LSA...9..56K [DOI: https://dx.doi.org/10.1038/s41377-020-0286-z] 56.

20. Maggiolini, E. et al. Strongly enhanced light–matter coupling of monolayer WS₂ from a bound state in the continuum. Nat. Mater. 22, 964–969 (2023).

21. Sortino, L. et al. Atomic-layer assembly of ultrathin optical cavities in van der Waals heterostructure metasurfaces. Nat. Photon (2025). https://doi.org/10.1038/s41566-025-01675-4.

22. Sortino, L et al. Enhanced light-matter interaction in an atomically thin semiconductor coupled with dielectric nano-antennas. Nat. Commun.; 2019; 10, 2019NatCo.10.5119S [DOI: https://dx.doi.org/10.1038/s41467-019-12963-3] 5119.

23. Benimetskiy, F et al. Measurement of local optomechanical properties of a direct bandgap 2D semiconductor. APL Mater.; 2019; 7, 101126.2019APLM..7j1126B [DOI: https://dx.doi.org/10.1063/1.5117259]

24. Kühner, L et al. Radial bound states in the continuum for polarization-invariant nanophotonics. Nat. Commun.; 2022; 13, 2022NatCo.13.4992K [DOI: https://dx.doi.org/10.1038/s41467-022-32697-z] 4992.

25. Mennel, L; Paur, M; Mueller, T. Second harmonic generation in strained transition metal dichalcogenide monolayers: MoS₂, MoSe₂, WS₂, and WSe₂. APL Photon.; 2019; 4, 034404.2019APLP..4c4404M [DOI: https://dx.doi.org/10.1063/1.5051965]

26. Lu, J et al. Periodically poled thin-film lithium niobate microring resonators with a second-harmonic generation efficiency of 250,000%/W. Optica; 2019; 6, pp. 1455-1460.2019Optic..6.1455L [DOI: https://dx.doi.org/10.1364/OPTICA.6.001455]

27. Fedotova, A et al. Second-harmonic generation in resonant nonlinear metasurfaces based on lithium niobate. Nano Lett.; 2020; 20, pp. 8608-8614.2020NanoL.20.8608F [DOI: https://dx.doi.org/10.1021/acs.nanolett.0c03290]

28. Karvounis, A; Timpu, F; Vogler-Neuling, VV; Savo, R; Grange, R. Barium titanate nanostructures and thin films for photonics. Adv. Opt. Mater.; 2020; 8, 2001249. [DOI: https://dx.doi.org/10.1002/adom.202001249]

29. Xu, L et al. Forward and backward switching of nonlinear unidirectional emission from GaAs nanoantennas. ACS Nano; 2019; 14, pp. 1379-1389. [DOI: https://dx.doi.org/10.1021/acsnano.9b07117]

30. Tilmann, B et al. Comparison of harmonic generation from crystalline and amorphous gallium phosphide nanofilms. Adv. Opt. Mater.; 2023; 11, 2300269. [DOI: https://dx.doi.org/10.1002/adom.202300269]

31. Smirnova, D; Kivshar, YS. Multipolar nonlinear nanophotonics. Optica; 2016; 3, pp. 1241-1255.2016Optic..3.1241S [DOI: https://dx.doi.org/10.1364/OPTICA.3.001241]

32. Kivshar, Y. All-dielectric meta-optics and non-linear nanophotonics. Natl. Sci. Rev.; 2018; 5, pp. 144-158. [DOI: https://dx.doi.org/10.1093/nsr/nwy017]

33. Sain, B; Meier, C; Zentgraf, T. Nonlinear optics in all-dielectric nanoantennas and metasurfaces: a review. Adv. Photon.; 2019; 1, 024002.2019AdPho..1b4002S [DOI: https://dx.doi.org/10.1117/1.AP.1.2.024002]

34. Löchner, FJ et al. Polarization-dependent second harmonic diffraction from resonant GaAs metasurfaces. ACS Photonics; 2018; 5, pp. 1786-1793. [DOI: https://dx.doi.org/10.1021/acsphotonics.7b01533]

35. Carletti, L et al. Steering and encoding the polarization of the second harmonic in the visible with a monolithic LiNbO₃ metasurface. ACS Photonics; 2021; 8, pp. 731-737. [DOI: https://dx.doi.org/10.1021/acsphotonics.1c00026]

36. Wang, C et al. Metasurface-assisted phase-matching-free second harmonic generation in lithium niobate waveguides. Nat. Commun.; 2017; 8, 2017NatCo..8.2098W [DOI: https://dx.doi.org/10.1038/s41467-017-02189-6] 2098.

37. Kuo, PS; Bravo-Abad, J; Solomon, GS. Second-harmonic generation using-quasi-phasematching in a GaAs whispering-gallery-mode microcavity. Nat. Commun.; 2014; 5, 2014NatCo..5.3109K [DOI: https://dx.doi.org/10.1038/ncomms4109] 3109.

38. Santiago-Cruz, T et al. Photon pairs from resonant metasurfaces. Nano Lett.; 2021; 21, pp. 4423-4429.2021NanoL.21.4423S [DOI: https://dx.doi.org/10.1021/acs.nanolett.1c01125]

39. Hsu, CW et al. Observation of trapped light within the radiation continuum. Nature; 2013; 499, pp. 188-191.2013Natur.499.188H [DOI: https://dx.doi.org/10.1038/nature12289]

40. Plotnik, Y et al. Experimental observation of optical bound states in the continuum. Phys. Rev. Lett.; 2011; 107, 183901.2011PhRvL.107r3901P [DOI: https://dx.doi.org/10.1103/PhysRevLett.107.183901]

41. Rybin, MV et al. High-Q supercavity modes in subwavelength dielectric resonators. Phys. Rev. Lett.; 2017; 119, 243901.2017PhRvL.119x3901R [DOI: https://dx.doi.org/10.1103/PhysRevLett.119.243901]

42. Hsu, CW; Zhen, B; Stone, AD; Joannopoulos, JD; Soljačić, M. Bound states in the continuum. Nat. Rev. Mater.; 2016; 1, pp. 1-13. [DOI: https://dx.doi.org/10.1038/natrevmats.2016.48]

43. Koshelev, K; Bogdanov, A; Kivshar, Y. Meta-optics and bound states in the continuum. Sci. Bull.; 2019; 64, pp. 836-842. [DOI: https://dx.doi.org/10.1016/j.scib.2018.12.003]

44. Koshelev, K; Lepeshov, S; Liu, M; Bogdanov, A; Kivshar, Y. Asymmetric metasurfaces with high-Q resonances governed by bound states in the continuum. Phys. Rev. Lett.; 2018; 121, 193903.2018PhRvL.121s3903K [DOI: https://dx.doi.org/10.1103/PhysRevLett.121.193903]

45. Carletti, L; Koshelev, K; De Angelis, C; Kivshar, Y. Giant nonlinear response at the nanoscale driven by bound states in the continuum. Phys. Rev. Lett.; 2018; 121, 033903.2018PhRvL.121c3903C [DOI: https://dx.doi.org/10.1103/PhysRevLett.121.033903]

46. Liu, Z et al. High-Q quasibound states in the continuum for nonlinear metasurfaces. Phys. Rev. Lett.; 2019; 123, 253901.2019PhRvL.123y3901L [DOI: https://dx.doi.org/10.1103/PhysRevLett.123.253901]

47. Koshelev, K et al. Subwavelength dielectric resonators for nonlinear nanophotonics. Science; 2020; 367, pp. 288-292.2020Sci..367.288K [DOI: https://dx.doi.org/10.1126/science.aaz3985]

48. Parry, M et al. Enhanced generation of nondegenerate photon pairs in nonlinear metasurfaces. Adv. Photon.; 2021; 3, 055001.2021AdPho..3e5001P [DOI: https://dx.doi.org/10.1117/1.AP.3.5.055001]

49. Santiago-Cruz, T et al. Resonant metasurfaces for generating complex quantum states. Science; 2022; 377, pp. 991-995.2022Sci..377.991S [DOI: https://dx.doi.org/10.1126/science.abq8684]

50. Sharapova, PR; Kruk, SS; Solntsev, AS. Nonlinear dielectric nanoresonators and metasurfaces: Toward efficient generation of entangled photons. Laser Photon. Rev.; 2023; 17, 2200408.2023LPRv..1700408S [DOI: https://dx.doi.org/10.1002/lpor.202200408]

51. Zograf, G et al. High-harmonic generation from resonant dielectric metasurfaces empowered by bound states in the continuum. ACS Photonics; 2022; 9, pp. 567-574. [DOI: https://dx.doi.org/10.1021/acsphotonics.1c01511]

52. Zalogina, A et al. High-harmonic generation from a subwavelength dielectric resonator. Sci. Adv.; 2023; 9, eadg2655. [DOI: https://dx.doi.org/10.1126/sciadv.adg2655]

53. Green, TD et al. Optical material anisotropy in high-index transition metal dichalcogenide Mie nanoresonators. Optica; 2020; 7, pp. 680-686.2020Optic..7.680G [DOI: https://dx.doi.org/10.1364/OPTICA.389192]

54. Munkhbat, B; Wróbel, P; Antosiewicz, TJ; Shegai, TO. Optical constants of several multilayer transition metal dichalcogenides measured by spectroscopic ellipsometry in the 300–1700 nm range: high index, anisotropy, and hyperbolicity. ACS Photonics; 2022; 9, pp. 2398-2407. [DOI: https://dx.doi.org/10.1021/acsphotonics.2c00433]

55. Verre, R et al. Transition metal dichalcogenide nanodisks as high-index dielectric mie nanoresonators. Nat. Nanotechnol.; 2019; 14, pp. 679-683.2019NatNa.14.679V [DOI: https://dx.doi.org/10.1038/s41565-019-0442-x]

56. Munkhbat, B; Küçüköz, B; Baranov, DG; Antosiewicz, TJ; Shegai, TO. Nanostructured transition metal dichalcogenide multilayers for advanced nanophotonics. Laser Photon. Rev.; 2023; 17, 2200057.2023LPRv..1700057M [DOI: https://dx.doi.org/10.1002/lpor.202200057]

57. Weber, T et al. Intrinsic strong light-matter coupling with self-hybridized bound states in the continuum in van der Waals metasurfaces. Nat. Mater.; 2023; 22, pp. 970-976.2023NatMa.22.970W [DOI: https://dx.doi.org/10.1038/s41563-023-01580-7]

58. Bouteyre, P et al. Simultaneous observation of bright and dark polariton states in subwavelength gratings made from quasi-bulk WS₂. Opt. Express; 2025; 33, pp. 21173-21182. [DOI: https://dx.doi.org/10.1364/OE.553437]

59. Busschaert, S et al. Transition metal dichalcogenide resonators for second harmonic signal enhancement. ACS Photonics; 2020; 7, pp. 2482-2488. [DOI: https://dx.doi.org/10.1021/acsphotonics.0c00751]

60. Panmai, M; Xiang, J; Zhou, L; Li, S; Lan, S. Revealing Mie resonances with enhanced and suppressed second-order nonlinear optical responses in a hexagonal-prism-like MoS₂ nanoparticle. Laser Photon. Rev.; 2023; 17, 2300346.2023LPRv..1700346P [DOI: https://dx.doi.org/10.1002/lpor.202300346]

61. Sung, J et al. Room-temperature continuous-wave indirect-bandgap transition lasing in an ultra-thin WS₂ disk. Nat. Photon.; 2022; 16, pp. 792-797.2022NaPho.16.792S [DOI: https://dx.doi.org/10.1038/s41566-022-01085-w]

62. Gan, X-T et al. Microwatts continuous-wave pumped second harmonic generation in few-and mono-layer GaSe. Light Sci. Appl.; 2018; 7, pp. 17126-17126. [DOI: https://dx.doi.org/10.1038/lsa.2017.126]

63. Zhao, M et al. Atomically phase-matched second-harmonic generation in a 2D crystal. Light Sci. Appl.; 2016; 5, pp. e16131-e16131. [DOI: https://dx.doi.org/10.1038/lsa.2016.131]

64. Shi, J et al. 3R MoS₂ with broken inversion symmetry: a promising ultrathin nonlinear optical device. Adv. Mater.; 2017; 29, 1701486. [DOI: https://dx.doi.org/10.1002/adma.201701486]

65. Xu, X et al. Towards compact phase-matched and waveguided nonlinear optics in atomically layered semiconductors. Nat. Photon.; 2022; 16, pp. 698-706.2022NaPho.16.698X [DOI: https://dx.doi.org/10.1038/s41566-022-01053-4]

66. Hallil, H et al. Strong piezoelectricity in 3R-MoS₂ flakes. Adv. Electron. Mater.; 2022; 8, 2101131. [DOI: https://dx.doi.org/10.1002/aelm.202101131]

67. Dong, Y et al. Giant bulk piezophotovoltaic effect in 3R-MoS₂. Nat. Nanotechnol.; 2023; 18, pp. 36-41.2023NatNa.18..36D [DOI: https://dx.doi.org/10.1038/s41565-022-01252-8]

68. Zograf, G et al. Combining ultrahigh index with exceptional nonlinearity in resonant transition metal dichalcogenide nanodisks. Nat. Photon.; 2024; 18, pp. 751-757. [DOI: https://dx.doi.org/10.1038/s41566-024-01444-9]

69. Weissflog, MA et al. A tunable transition metal dichalcogenide entangled photon-pair source. Nat. Commun.; 2024; 15, [DOI: https://dx.doi.org/10.1038/s41467-024-51843-3] 7600.

70. Qin, B et al. Interfacial epitaxy of multilayer rhombohedral transition-metal dichalcogenide single crystals. Science; 2024; 385, pp. 99-104. [DOI: https://dx.doi.org/10.1126/science.ado6038]

71. Feng, J et al. Polarization-entangled photon-pair source with van der Waals 3R-WS₂ crystal. eLight; 2024; 4, 16. [DOI: https://dx.doi.org/10.1186/s43593-024-00074-6]

72. Trovatello, C et al. Quasi-phase-matched up-and down-conversion in periodically poled layered semiconductors. Nat. Photon.; 2025; 19, pp. 291-299. [DOI: https://dx.doi.org/10.1038/s41566-024-01602-z]

73. Ling, H et al. Nonlinear van der Waals metasurfaces with resonantly enhanced light generation. Nano Lett.; 2025; 25, ASAP. [DOI: https://dx.doi.org/10.1021/acs.nanolett.5c00952]

74. Munkhbat, B et al. Transition metal dichalcogenide metamaterials with atomic precision. Nat. Commun.; 2020; 11, 2020NatCo.11.4604M [DOI: https://dx.doi.org/10.1038/s41467-020-18428-2] 4604.

75. Kühne, J et al. Fabrication robustness in BIC metasurfaces. Nanophotonics; 2021; 10, pp. 4305-4312. [DOI: https://dx.doi.org/10.1515/nanoph-2021-0391]

76. Liu, Q et al. Boosting second harmonic generation by merging bound states in the continuum. Phys. Rev. B; 2023; 107, 245408.2023PhRvB.107x5408L [DOI: https://dx.doi.org/10.1103/PhysRevB.107.245408]

77. Anthur, AP et al. Continuous wave second harmonic generation enabled by quasi-bound-states in the continuum on gallium phosphide metasurfaces. Nano Lett.; 2020; 20, pp. 8745-8751.2020NanoL.20.8745A [DOI: https://dx.doi.org/10.1021/acs.nanolett.0c03601]

78. Lee, J et al. Observation and differentiation of unique high-Q optical resonances near zero wave vector in macroscopic photonic crystal slabs. Phys. Rev. Lett.; 2012; 109, 067401.2012PhRvL.109f7401L [DOI: https://dx.doi.org/10.1103/PhysRevLett.109.067401]

79. Dewambrechies, A; Polyakov, AY; Küçüköz, B; Shegai, TO. Enhanced second-order nonlinearities at strained ultrasharp zigzag edges in multilayer MoS₂. J. Phys. Chem. C; 2023; 127, pp. 15395-15405. [DOI: https://dx.doi.org/10.1021/acs.jpcc.3c03812]

80. Polyakov, AY et al. Top-down nanostructured multilayer MoS₂ with atomically sharp edges for electrochemical hydrogen evolution reaction. Mater. Today Nano; 2024; 25, 100467. [DOI: https://dx.doi.org/10.1016/j.mtnano.2024.100467]

81. Koshelev, K; Toftul, I; Hwang, Y; Kivshar, Y. Scattering matrix for chiral harmonic generation and frequency mixing in nonlinear metasurfaces. J. Opt.; 2024; 26, 055003.2024JOpt..26e5003K [DOI: https://dx.doi.org/10.1088/2040-8986/ad3a78]

82. Lobato, I; Friedrich, T; Van Aert, S. Deep convolutional neural networks to restore single-shot electron microscopy images. npj Comput. Mater.; 2024; 10, 10.2024npjCM.10..10L [DOI: https://dx.doi.org/10.1038/s41524-023-01188-0]

83. O’Brien, K et al. Predicting nonlinear properties of metamaterials from the linear response. Nat. Mater.; 2015; 14, pp. 379-383.2015NatMa.14.379O [DOI: https://dx.doi.org/10.1038/nmat4214]

84. Zu, R et al. Analytical and numerical modeling of optical second harmonic generation in anisotropic crystals using #SHAARP package. npj Comput. Mater.; 2022; 8, 246.2022npjCM..8.246Z [DOI: https://dx.doi.org/10.1038/s41524-022-00930-4]

85. Psilodimitrakopoulos, S; Ilin, S; Zelenkov, LE; Makarov, S; Stratakis, E. Tailoring of the polarization-resolved second harmonic generation in two-dimensional semiconductors. Nanophotonics; 2024; 13, pp. 3181-3206. [DOI: https://dx.doi.org/10.1515/nanoph-2024-0267]

86. Zograf, G. & Antosiewicz, T. J. Data and code for “Ultrathin 3R-MoS₂ metasurfaces with atomically precise edges for efficient nonlinear nanophotonics”. Zenodohttps://doi.org/10.5281/zenodo.14194107 (2025).