1.

Introduction

Metasurfaces are thin and flat arrays of subwavelength nanoparticles, enabling control over the polarization, phase, amplitude, and dispersion of light.¹ They can be used for light emission, detection, modulation, control, and/or amplification at the nanoscale. In recent years, metasurfaces have been a subject of ongoing intense studies as their optical properties can be adapted to a diverse set of applications, including superlenses, tunable images, holograms, etc.^1–5 High-refractive-index dielectric metasurfaces provide a powerful platform for controlling light that can go beyond plasmonics as they cause negligible losses as compared with plasmonic metasurfaces. Dielectric materials offer the unique ability to efficiently manipulate light at the nanoscale based on the simultaneous excitation and control over the optically induced electric and magnetic Mie-type resonances. Resonant dielectric metasurfaces with high-quality factor (high-$Q$) resonances, in particular, are of significant interest due to their ability to strongly enhance the electromagnetic near-fields and boost the light–matter interactions at the nanoscale. In other words, they allow for enhancing the response of metasurfaces to an external electromagnetic field at a particular frequency. Moreover, high-$Q$ metasurfaces can increase the storage time of photons and, thus, light–matter interactions within the subwavelength resonators. This will facilitate various nanophotonics applications, such as enhanced nonlinear photon generation, optical sensing, optoacoustic vibrations, and narrowband filtering. In the last decade, high-$Q$ metasurfaces were mainly associated with Fano resonances (FRs) featuring asymmetric spectral line profiles.^6–10 In FRs, the asymmetry originates from a close interaction of a discrete (localized) state with a continuum of propagation modes.^11–13 Some examples are trapped modes in arrays of dielectric nanodisks with asymmetric holes,^14–16 ring and disk cavities,¹⁷ Dolmen structures,^18,19 and aggregated nanoparticles.^20–23 Recently, different approaches for generating high-$Q$ resonances have been based on the bound-state-in-the-continuum (BIC), i.e., a localized state with zero linewidth that is embedded in the continuum.^15,24–32 Indeed, optical BICs provide a unique opportunity to manipulate the light–matter interaction within the radiative continuum because of their ultrahigh-$Q$ origin and associated giant enhancement of the electromagnetic near-field.³³

By contrast, designing metasurfaces with high-$Q$ resonances is usually achieved via continuous tuning of parameters, with limited control on the linewidth, amplitude, and spectral positions. Currently, one of the typical approaches for designing metasurfaces with the desired resonance is based on direct optimization of one or two parameters via brute-force simulations.^34–44 This is a time-consuming task accompanied by a random success on the output parameters of the desired resonances.⁴⁵ Recently, deep-learning approaches, based on artificial neural networks (ANNs), have emerged as a revolutionary and robust methodology in nanophotonics.^46–60 Indeed, applying deep learning algorithms to the nanophotonic inverse design can introduce remarkable design flexibility that can exceed that of the conventional methods.⁶¹ The inverse-design approach is based on a training process that enables fast prediction of complex optical properties of nanostructures with intricate architectures.

As a (nonunique) example, we have targeted toroidal dipoles (TDs) due to their promising applications in the formation of anapole states and electromagnetic energy localization. We investigate the nonradiating TD supported by two parallel silicon bars as the building blocks of the metasurface. The reason for choosing such a geometry is that it allows for a reasonable number of parameters to be optimized as a proof of concept. The parameters include the length and width of the bars, as well as the gap between them. At the same time, this TD corresponds to the symmetry-protected BIC. Subsequently, with this TD-BIC model, we demonstrate the deep-learning-assisted inverse design of arbitrary high-$Q$ resonances with different line widths, amplitudes, and spectral positions. We employ a multilayer perceptron (MLP) variant of the ANN as our model.^62–64 An MLP consists of multiple layers of perceptrons, including an input layer and an output layer with several hidden layers. In this work, each artificial neuron in one layer connects with a certain weight to every neuron in the following layers that are adapted during the learning state. Once learned, the weight values remain fixed, and the model can be used to infer the target metasurface design parameters. The hidden layers establish a nonlinear mapping between the input and output through training from the given dataset, which is then able to predict the response of the system or inversely determine the design parameters for the desired performance.

The proposed deep-learning-assisted inverse-design approach provides a platform to design metasurfaces for more than one application. In this paper, we employ machine learning to obtain a bifunctional metasurface dealing with photons and phonons, simultaneously.^65–74 Photon–photon conversions, so-called nonlinear nanophotonics, is at the heart of modern macroscopic optics, including lasers, sensors, imaging, and information technology. By contrast, photon–phonon conversions are the state-of-the-art solution for precision mass sensors, micromanipulation, and sensing biochemical materials, with transformative implications in the fields of health and security. A combined photon–phonon conversion can be used for nonionizing and noninvasive imaging.⁷³ Here, using the deep-learning-assisted inverse approach, we design and fabricate a single optoacoustic metasurface that enhances third-harmonic generation (THG) intensity up to 400 times and acoustic mode excitations up to 100 times, concurrently, all through a designed high-$Q$ resonance associated with a strong electric near-field enhancement. The inverse-design approach proposed in this paper is extendable to other characteristics and applications of metasurfaces and significantly circumvents the time-consuming, case-by-case numerical models in conventional electromagnetic nanostructure designs.

2.

Results and Discussions

To obtain the initial high-$Q$ resonances, we have defined the building blocks of metasurfaces to be two identical silicon nanobars with width $w$, length $L$, and the offset $x_0$, which is the distance between the center of the two bars fabricated on a glass substrate, as shown in the top panel of Fig. 1(a). It is worth mentioning that a large variety of other geometries, for generating high-$Q$ resonances have been demonstrated.^{15,24–31,35,39} However, discussions on advantages and disadvantages of various geometries are beyond the scope of the current study. Instead, we concentrate on customizing the generated high-$Q$ resonances. The thickness of bars is fixed at 150 nm, and the periodicity of unit cells is fixed at $D=900\mathrm{nm}$ in both the $x$ and $y$ directions. As can be seen in the bottom panel of Fig. 1(a), this structure supports a strong FR in the transmission spectrum when the incident light is polarized along the $y$ axis. This FR is formed by the interference and coupling between a “bright” electric dipole (ED) resonance $p_y$ and a “dark” TD mode $T_y$. Figure 1(b) shows the corresponding spherical multipolar decomposition of the metasurface. As can be seen, the optical response is dominated by ED excitation with a small contribution from the magnetic quadrupole (MQ) resonance. Such a pronounced ED feature was further investigated by performing the Cartesian multipolar analysis [Fig. 1(c)]. It is worth noting that the ED response is mainly due to the strong excitation of TD mode $T_y$ with an in-plane ED mode $p_y$, which is polarized along the same direction as the optical pump.

Fig. 1

(a) Top: (top left) Schematics of the silicon nanobars metasurface and (top right) its unit cell. Bottom: Calculated transmission spectrum of the metasurface with structural parameters $w=316\mathrm{nm}$, $L=580\mathrm{nm}$, $x_0=189\mathrm{nm}$. (b) Spherical multipolar structure of the metasurface. (c) Top: Cartesian ED and TD modes excitations. Bottom: The electric energy enhancement ${\eta }_E/{\eta }_{E_0}$. It is defined as the electric energy inside the two nanobars normalized by the electric energy within the same volume of the nanobars for the pump field. (d) Electric near-field distributions at the resonance. Left: 3-D view. Right: top view.

Interestingly, due to the $C_2$ symmetry of our subdiffractive system, the TD and MQ do not contribute to the far-field radiation along the $z$ direction. The far-field optical response is dominated by the ED mode $p_y$. The nonradiating TD mode is a symmetry-protected BIC, where the ED mode $p_y$ plays a role in opening a leaky channel and transforming this ideal BIC into a quasi-BIC with a finite $Q$-factor. In our previous work,¹⁵ we studied the formation of BIC enabled by a magnetic dipole (MD) resonance in Si disk-with-hole metasurfaces, where geometrical asymmetry was introduced to open a leaky channel. In contrast, here we show that the leaky channel, i.e., the excitation of ED mode $p_y$, can be formed directly by properly choosing the structural dimensions of the symmetric nanobars. Owing to the nonradiating nature of the dominant resonance—TD mode, a clear enhancement of the stored electric energy inside the nanobars is observed, as shown in the bottom of Fig. 1(c). Figure 1(d) gives the calculated electric near-field distributions. A pronounced poloidal current distribution can be observed from the two nanobars, indicating strong TD excitation. The small portion of MQ excitation shown in Fig. 1(d) is due to the uncompensated circulating magnetic field in the flat geometry, formed by two antiparallel MD moments at the nodes of the poloidal current distribution. A comparison between the electric near-fields between the $y$-polarized pump and $x$-polarized pump incidence can be seen in Fig. S1 of the Supplemental Material. We observed significant near-field enhancement inside the nanostructure for the $y$-polarized pump incidence as compared with the case for the $x$-polarized pump incidence. A further investigation on the band structure and the corresponding mode profiles can be seen in Figs. S2 and S3 of Section I in the Supplemental Material.

As mentioned earlier, achieving scalable metasurfaces with several interdependent characteristics, including quality factor and spectral position, is the main target of this paper. There is a significant demand in the photonics community to achieve both conditions rather than solely obtaining a high-$Q$ resonance. In this respect, we employ the deep-learning approach to inversely design our metasurfaces while simultaneously controlling and optimizing on $Q$-factor, amplitude, and spectral position. For this task, we use the open-source neural-network library Keras⁷⁵ written in Python to implement our method.

For the inverse design of nanophotonic structures using a deep-learning approach, one challenge is that the same far-field electromagnetic response can correspond to different designs, i.e., several different structures can give the same responses. This nonuniqueness of the response-to-design mapping will induce conflicting examples within the training set and may prevent convergence. To avoid this issue, we use the tandem network (TN) approach,^46,56 as shown in Fig. 2. The TN architecture consists of an inverse-design network connected to a forward model network. The forward network learns the mapping from the structural parameters to the optical responses and is trained separately first. After the forward network is trained, it is placed after the inverse-design model network, and its network weights remain fixed during the training of the inverse-design model network. The inverse-design network learns a mapping from the optical responses to the structural parameters. When training the inverse-design network, its weights are updated to minimize the loss objective: $J=\frac{1}{n}\sum _i{(r_i-o_i)}^2$ with $r_i$ and $o_i$ being the values predicted by the neural network and the ground truth of the response (forward model network) or the structural parameters (inverse-design model network). Trained in this way, the inverse-design network is not constrained to produce a prespecified design. Instead, it is free to infer any design that results in the desired forward behavior.

Fig. 2

The architecture of the TN model, which consists of an inverse-design network connected to a pretrained forward model network. $X$ represents the input and output, which is the transmission spectra data in our case, and $Y$ represents the output in the middle layer which is the structural parameters here.

In our case, for simplicity, we specify three structural parameters of the nanobars to learn: width $w$, length $L$, and the offset distance $x_0$. It is worth mentioning that the machine-learning approach can be extended to predict nanostructures with more parameters or even different types of parameters such as material properties and material losses. By randomly specifying them, we first use the rigorous coupled-wave analysis (RCWA)⁷⁶ to generate 25,000 training examples, in which we obtain the transmission spectra of the metasurfaces covering a wavelength range of 1400 to 1600 nm. Here, the refractive index of silicon in this wavelength range is nearly constant; thus in our simulation, we keep the refractive index of silicon as 3.6 to expedite the training data generation process. RCWA is a frequency-domain modal method based on the decomposition of the periodic structure and the pseudoperiodic solution of Maxwell’s equations in terms of their Fourier expansions.⁷⁷ It has been widely used for modeling light responses from periodic optical structures due to its fast convergence and accurate far-field calculations.⁷⁸ It is quite suitable for modeling the electromagnetic responses of metasurfaces and generating massive training data, especially when considering the inverse design based on neural networks.^46,51 It is worth mentioning that, for multilayer structures, the transfer matrix method can be one more choice.⁷⁹

The forward model network is designed to have four fully connected layers with each layer having 400–600–400–200 dimensions, respectively. We set the learning parameter batch size as 256, and use a learning rate of 0.001 and decay of $1\times {10}^{-6}$. We first train the forward model network and evaluate it to see how well it can predict the given transmission spectra. Figure 3(a) shows the learning curves for training and validation loss. It can be seen that both the training loss and validation loss decrease significantly after 10,000 epochs of training and become $<0.005$ after 30,000 epochs of training. This indicates that the trained network can estimate an appropriate spectrum that is similar to the spectrum calculated analytically. As a test example, we use RCWA to simulate the transmission spectrum for an individual Si nanobar metasurface with parameters $[w,L,x_0]$ being [300, 700, 300] nm, as shown by the black dashed curve of Fig. 3(b). Then, we input these structural parameters to the network and predict the corresponding output, which is shown in the blue curve of Fig. 3(b). As can be seen, our forward network can predict the transmission spectrum from our metasurface accurately. We have also tested our forward network by specifying different structural parameter sets corresponding to different transmission spectra (see Fig. S4 in the Supplemental Material), which all verify the effectiveness of our forward network.

Fig. 3

(a) Evolution of the training loss for the forward model network. (b) Comparison of the NN approximation to the real transmission spectrum. (c) Evolution of the training loss for the inverse-design model network. (d) Comparison of the spectra between the NN approximation and the input based on Eq. (2).

Next, we train the inverse-design model network by fixing the weights in the pretrained forward model network. Since the forward model network is differentiable, we can train the inverse-design model network with a loss placed after the forward model network. As mentioned above, this will overcome the issue of nonuniqueness in the inverse spectrum of electromagnetic waves, as the design by the neural network is not required to be identical to the design parameters that produced the training samples, but the spectrum inferred by the forward model network must match the target spectrum. The loss can be further lowered when the generated design and the real design have similar responses after training. The inverse-design model network has five fully connected layers with 600, 600, 400, 200, 200 dimensions, respectively. For the inverse design of our metasurfaces, the transmission spectrum is considered the input of the TN. The design parameters are predicted from the intermediate layer of the whole network. The training process is shown in Fig. 3(c). As can be seen, using the TN approach, the learning of inverse design has converged effectively. We then test our inverse-design network using a Fano formula to define the transmission spectrum:^80,81

We first design a target FR at ${\lambda }_0=1500\mathrm{nm}$, with resonance linewidth $\mathrm{\Delta }\lambda =5\mathrm{nm}$ and $q=0.5$, as shown by the black dashed curve of Fig. 3(d). Here, we input the resonance frequency linewidth $\mathrm{\Gamma }$ by the defined resonance wavelength linewidth $\mathrm{\Delta }\lambda$ by $\mathrm{\Gamma }=\mathrm{\Delta }\lambda {\omega }_0/{\lambda }_0$. We then use the network to predict the structural parameters $[w,L,x_0]$ of the required metasurface as $[\mathrm{316,580,189}]\mathrm{nm}$. The transmission spectrum of the predicted metasurface is shown in the dashed red curve of Fig. 3(d). It matches well with the desired Fano-shape curve based on Eqs. (1) and (2).

The deep network provides a powerful approach to design nanophotonic structures inversely. Once the training process is finished, the inverse design calculation takes only around 0.05 s in our case [only using CPU in a normal desktop with 64-bit Operating System: Processor: Intel(R) Core(TM) i7-4770 CPU @ 3.40 GHz, RAM: 16.0 GB], which is both promising and much more effective, as compared with conventional electromagnetic solvers. For nanophotonics applications, as mentioned above, it is a requirement to obtain scalable metasurfaces with controllable characteristics. In the following, using the TN approach, we further test the inverse design of FRs with different spectral positions, linewidths, or amplitudes. According to Eqs. (1) and (2), we first specify target FRs at ${\lambda }_0=1450$, 1500, and 1550 nm, respectively. Subsequently, we keep the linewidth $\mathrm{\Delta }\lambda =15\mathrm{nm}$ and $q=0.8$, as shown in the black dashed curve of Figs. 4(a)–4(c). We then use the trained neural network to predict the structural parameters of the metasurfaces that can provide such spectra. The transmission spectra of the designed metasurfaces are shown in the red curves of Figs. 4(a)–4(c). They match and satisfy the design goal well. By varying the value of parameter $q$ or the linewidth of the resonance $\mathrm{\Delta }\lambda$, the neural network can easily predict the metasurface design for the required FR with different characteristics, as shown in Figs. 4(d)–4(i). It provides a powerful method to control the near-field and electric energy confinement at the nanoscale (see Fig. S5 in the Supplemental Material).

Fig. 4

Inverse design of Si nanobar metasurfaces with Fano-shape transmission spectra. (a)–(c) ${\lambda }_0=1450$, 1500, and 1550 nm, respectively. $\mathrm{\Delta }\mathit{\lambda }=15\mathrm{nm}$, $q=0.8$. (d)–(f) ${\lambda }_0=1500\mathrm{nm}$, $\mathrm{\Delta }\mathit{\lambda }=10\mathrm{nm}$, $q=0.3$, 0.5, and 0.7, respectively. (g)–(i) ${\lambda }_0=1500\mathrm{nm}$, $\mathrm{\Delta }\mathit{\lambda }=5$, 15, and 25 nm, respectively, $q=0.7$.

So far, based on the deep-learning approach, we have optimized the design to obtain high-$Q$ resonances via altering various parameters, including the linewidth, spectral position, and amplitude, simultaneously. Such a high-$Q$ resonance can significantly facilitate nanostructures to enhance light–matter interactions for various applications such as nonlinear optics and optomechanics. As an example, here, we investigate the THG from three designed metasurfaces with resonances at 1450, 1500, and 1550 nm, as shown in Figs. 4(a)–4(c), respectively. Here, we have taken into account that the spectral full-width-at-half-maximum of our experimentally used laser is around 15 nm for the wavelength range 1400 to 1600 nm. Thus, we consider the resonances with this linewidth to maximize the nonlinear signal generation.

The fabrication is done via the standard electrolithography, similar to our previous work,¹⁵ in which nanostructures are fabricated out of amorphous silicon on a glass substrate. The scanning electron microscope (SEM) image of one fabricated metasurface sample with a designed resonance at 1500 nm is shown in Fig. 5(a). The measured dimensions of the fabricated metasurfaces $[w,L,x_0]$ are $[\mathrm{269,805,188}]\mathrm{nm}$ (resonant at 1450 nm), $[\mathrm{290,765,184}]\mathrm{nm}$ (resonant at 1500 nm), and $[\mathrm{310,780,184}]\mathrm{nm}$ (resonant at 1550 nm) (see Fig. S6 in the Supplemental Material). While the experimental and theoretical $x_0$ are identical, controlled by the lithography software, there are slight differences in the values of the fabricated $w$ and $L$ from the theoretical targeting values, as shown in Figs. 4(a)–4(c), due to fabrication imperfections. Such discrepancies, together with a possible minor inhomogeneity of the experimentally deposited silicon film, slightly affect the experimentally fabricated samples. It is worth noting that the lattice perturbations at the array’s edge can also break the periodic boundary conditions leading to scattering of light in all directions into the free space, which weakens the total transmission or total reflection of the FR.

Fig. 5

(a) SEM image of the fabricated sample with designed resonance at 1500 nm. (b) Experimentally measured linear spectra. (c) Experimentally measured THG spectra of the samples.

We first measured the linear transmission spectra of the three metasurfaces under plane wave normal incidence with the electric field polarized along the $y$ axis, as shown in Fig. 5(b). Pronounced asymmetric FRs were observed around the desired spectral positions. Subsequently, we performed the TH spectroscopy measurement. A femtosecond laser beam with 200 fs pulse width and 80 MHz repetition rate was focused by an aspheric lens with a focal length of 5 cm to a beam waist of $20\mu \mathrm{m}$. The pump polarization was adjusted along the $y$ axis to excite the designed TD BIC state and was tuned ranging from 1400 to 1600 nm, with maximum mean power in the sample plane up to around 66 mW, leading to a maximum peak intensity value of around ${0.66\mathrm{GW}/\mathrm{cm}}^2$. An objective with a numerical aperture (NA) $\mathrm{NA}=0.7$ was used to collect the transmitted TH emission power in the forward direction (see Fig. S7 in Section V. Experimental setup for nonlinear measurements of the Supplemental Material). The experimentally measured TH signals from the three designed metasurfaces are shown in Fig. 5(c). As can be seen, the TH signals are significantly enhanced around the resonances, while no THG enhancement is observed when the laser beam is polarized along the $x$ axis (see Fig. S8 in the Supplemental Material). By comparison, we observed 400-fold enhancement of the TH signal at the resonance position under the $y$-polarized pump as compared with the case of the $x$-polarized pump (Fig. S8 in the Supplemental Material). Similarly, by performing the nonlinear multipolar analysis, the TH signal is dominated by the TD excitation with small portions of MQ and EQ excitations, which exhibit the same $C_2$ symmetry (Fig. S9 in the Supplemental Material). This further leads to a stronger TH emission in the first-order diffraction compared with the zero-order diffraction due to the absence of coupling to these modes and the normal outgoing waves (see Fig. S9 in the Supplemental Material).

Another advantage of high-$Q$ resonances is that the particular electric field distributions can also facilitate the optomechanical vibration process.^82–86 In other words, by designing a high-$Q$ resonance through machine learning, one can engineer the vibrational modes of nanoparticles based on the optomechanical processes, too. Here, by considering the designed metasurface shown in Fig. 2(d), we numerically investigate the optomechanical vibrations and mechanical mode excitations. The mechanical vibration is modeled as

Due to the symmetric distributions of optical fields in these two nanobars shown in Figs. 1(d) and 6(a), the corresponding vibration process for them is similar and has the same absolute values. Subsequently, we analyze the vibrations by measuring the changes in the width and length along the center of right nanobar in the unit cell [see Fig. 6(a)]. The vibration displacements along the $x$ and $y$ directions, i.e., $D_x$ and $D_y$, respectively, are shown in Fig. 6(b). To reveal the excitation of acoustic modes under optical pump, we further calculate the spectral density of the optomechanical vibrations. The corresponding spectral density is computed by the Fourier transfer analysis of this time-dependent vibrational amplitude and is plotted in Fig. 6(c), which shows that the coherent phonon oscillation frequency for this nanobar is around 12 GHz for $D_x$ and 9.5 GHz for $D_y$. For comparison, Figs. 6(d)–6(f) show the corresponding optomechanical vibrations when the pump is polarized along the $x$ axis. The amplitude of optomechanical vibration is around 100 times stronger from the metasurface under $y$-polarized pump incidence, which corresponds to the excitation of the designed TD BIC state. It is worth noting that the optomechanical response and the excitation of acoustic modes are dependent on the dimensions of the nanostructures. Therefore, by geometric tuning of the nanostructures, efficient excitation of different acoustic modes at other frequencies can be achieved.

Fig. 6

(a)–(c) Optomechanic vibration under the $y$-polarized pump. (a) Displacement of the nanobars after 1 ns. (b) The transient displacement $D_x$ and $D_y$. (c) Spectral densities of displacement $D_x$ and $D_y$. (d)–(f) Optomechanical vibration under the $x$-polarized pump. (d) Displacement of the nanobars after 1 ns. (e) The transient displacement $D_x$ and $D_y$. (f) Spectral densities of displacement $D_x$ and $D_y$.

To demonstrate the importance of high-$Q$ resonance on acoustic modes, we further calculate the spectral densities for different pump wavelengths under $y$-polarized pump incidence, as shown in the three-dimensional (3-D) map in Fig. 7. As can be seen, a significant peak with the frequency of 12 GHz for $D_x$ (9.5 GHz for $D_y$) appears for the wavelength at 1500 nm due to the excitation of the TD BIC state at the optical pump. The excitation strength of the acoustic mode decreases dramatically when the optical pump is away from the TD BIC state. We then estimate the feedback of the mechanical vibration on the resonant optical response. Based on the transient vibration shown in Fig. 6(b), by assuming deformation around 50-pm displacement in the $x$ direction or 25-pm displacement in the $y$ direction, the maximum sensitivity of the scattering response in such metasurfaces is around 0.055% ${\mathrm{pm}}^{-1}$ (see Fig. S10 in the Supplemental Material). Utilizing the high-quality TD BIC state, up to 4.5% modulation of the transmission near the resonance can be expected when using a pulse laser with peak intensity $50\mathrm{GW}/{\mathrm{cm}}^2$ through the radiation force on the silicon nanostructures. These results suggest new opportunities for optomechanical applications such as light modulation and nanosensing with nanostructures.

Fig. 7

The spectral density of $D$ in the (a) $x$ and (b) $y$ directions for different laser pump wavelengths.

3.

Conclusion

Applying a machine-learning approach, we have demonstrated the inverse design of high-quality FR metasurfaces composed of two nanobars with scalable characteristics including the spectral position, line width, and amplitude of the transmission. The FR is originated from the TD-BIC state, featuring a strong near-field enhancement and intense electric energy localized inside the nanobars. We further employ these metasurfaces to simultaneously enhance photon–photon and photon–phonon interactions and achieve up to 400-fold THG enhancement and up to 100-fold enhancement of the amplitude of optomechanical vibrations. Our proposed scalable metasurfaces suggest new opportunities to control and enhance light–matter interactions, showing promising applications for realizing optoacoustic nonlinear metasurfaces.