Multipolar second-harmonic generation from high-Q quasi-BIC states in subwavelength resonators

2.1 Multipolar model of quasi-BIC formation

We consider a high-index cylindrical dielectric resonator that supports leaky modes (modes of an open resonator) that may hybridize (couple) when tuning geometric parameters. The particle is characterized by a frequency-dependent dielectric permittivity ε(ω) and is surrounded by a homogeneous host medium with εh=1. The harmonic time dependence of the fields in the form eiωt is implied.

Here, we focus on rotationally symmetric TE-polarized modes in the cylindrical coordinate system (Eφ,Hρ,Hz), whose electromagnetic field does not depend on the azimuthal angle φ. With the use of the finite element method (FEM) in COMSOL Multiphysics, we perform the eigenmode analysis numerically and plot dispersion as a function of the normalized geometric parameters defined as the disk aspect ratio r/h and the size parameter r/λ₀, where λ₀ is the mode wavelength. The size parameter defined as a ratio of the disk radius to the mode wavelength can also be referred to as a normalized eigenfrequency. The results of calculations for the AlGaAs nanodisk with the fixed height h = 645 nm are summarized in Figure 2. Two dispersion curves depicted as colored dots exhibit characteristic avoided crossing in the plane of parameters, which is a signature of the strong-coupling regime. It is accompanied by modification of the modes quality factors and formation of the quasi-BIC state that can be naturally understood through multipolar transformations of coupled modes as we describe below.

Figure 2:

Dispersion of the axially symmetric eigenmodes of the AlGaAs nanodisk. Colored circles depict the multipolar decomposition of these modes: magnetic dipole (MD – blue) and magnetic octupole (MO – pink), sizes of circles scale their relative contributions and radiative decay. The upper branch traces the high-quality mode and features a special point in the parameter space (r/h = 0.71, r/λ₀ = 0.29), where the quality factor of this mode reaches its maximum. The filled pink dot marks this quasi-BIC point, it corresponds to pure MO. Insets show near-field distributions of the electric field magnitude and far-field diagrams: the top row corresponds to the high-quality mode at r/h = 0.6, 0.71 (BIC), 0.8; the bottom row corresponds to the low-quality mode at r/h = 0.6, 0.71, 0.77. The corresponding values of aspect ratio are marked by dots at the horizontal axis. The overlaid black line sketches a dispersion branch of the mode near the second-harmonic frequency (the wavelength is doubled for direct comparison).

For a nonspherical shape of the resonator, the interacting modes can be in general viewed as a superposition of the spherical multipoles distinguished by orbital l and azimuthal m indices. Multipolar analysis of the considered TE-polarized axially symmetric eigenmodes suggests that the basis of parent multipoles mainly consists of a longitudinal magnetic dipole (MD) l = 1, m = 0 and a magnetic octupole (MO) l = 3, m = 0 so that each mode has two multipolar radiation channels, dipolar and octupolar. In Figure 2, sizes of circles illustrate their relative contributions in multipolar expansions of the modes (line plots are provided in Supplementary Material). The occurrence of high-Q supercavity mode is accompanied by an increase in the order of a dominating multipole from l = 1 (MD) to l = 3 (MO) and corresponds to the pure magnetic octupole. In the quasi-BIC condition (r/h = 0.71 and r/λ₀ = 0.29), two magnetic dipoles interfere destructively in the coupling to the octupole, thus, restoring its high-quality factor. Insets show how the electric field distributions inside the disk and far-field diagrams of these modes change as we move along the dispersion curves. While far from the BIC point, the radiation patterns are dipolar, near the BIC point the pattern of the high-Q mode turns to the rotationally symmetric magnetic octupole with three lobes. In addition, the black line in Figure 2 shows the dispersion of the mode near the second-harmonic frequency. The eigenmode analysis is performed for the AlGaAs disk taking into account the material dispersion. The refractive index in the fundamental wavelength range 1500–1700 nm n(ω)=10.73≈3.27 and in the SH wavelength range 750–850 nm n(2ω)≈3.52 [37].

We then corroborate our rigorous numerical results by an analytical model setting the correspondence of the non-Hermitian coupled-mode theory and the multipolar analysis. Coupling of two leaky eigenmodes shown in Figure 2 can be described by the Hamiltonian [33], [34]

being a square matrix, where E1=ω˜1+iγ1 and E2=ω˜2+iγ2 are complex frequencies of the modes 1 and 2, W and V are complex coupling parameters including the interaction between the modes via free space, W≠V*. A scheme of this phenomenological two-level model is shown in Figure 3A. The eigenfrequencies depend on the aspect ratio ξ≡r/h and can be deviated from the point ω˜1,2=0. Taking γ1=γ2=γ, V=W=ϰ+iγ, in the strong-coupling regime γ≲ϰ, the eigenvalues of the matrix 1 are given by E+=ϰ+2iγ and E−=−ϰ, that yields a frequency splitting of 2ϰ and the infinitely-large Q factor for one of the modes (BIC state). The avoided crossing of two modes in the two-level system qualitatively captures the numerically calculated dispersion depicted in Figure 2.

Figure 3:

(A) Schematic of interaction between two-coupled eigenmodes of an open dielectric resonator. Arrows illustrate individual couplings of the modes to the continuum and couplings between the modes via the continuum. (B) Schematic of the three-level multipolar model for the formation of the quasi-BIC state captured by the Hamiltonian Eq. (6).

For our system, we derive the dynamic equations for two modes with slowly varying amplitudes B1 and B2 excited by the external pump source in the following form:

where the matrix H^2 explicitly reads

The d coefficients govern both the dipolar and octupolar radiative decays of the modes and their coupling to each other and to the external field, d01,022≪d1,22. The near-field coupling is incorporated by off-diagonal ϰ terms. Using the scattering-matrix method (see Supplementary Material), the coupled-mode equations can be mapped onto the three-state model schematically illustrated in Figure 3B. The reduced multipolar basis includes two magnetic-dipole states a1,2 and the magnetic-octupole state a3 with m=0. The evolution equations recast as

with the amplitude column-vector a=[a1,a2,a3], and the external source

The non-Hermitian Hamiltonian of the structure in Eq. (4) assumes the form

where d1,2,0 are the effective dipole and octupole moments. The τ coefficients describe the hybridization of the dipole and octupole multipoles. Ohmic losses are neglected. The high-quality resonant state is formed at the frequency ωBIC=(ω1d2τ2+ω2d1τ1)/(d1τ1+d2τ2), where the imaginary part of the octupole-like eigenstate frequency is restored to a bare value d02.

The excitation of the quasi-BIC mode depends on how efficiently the incident source couples to the BIC state, that is on their spatial overlap. The external source can be divided into the octupolar component S0 and the dipolar component Sd and their fractions determine the resultant response. We consider two types of pump radiation which carry Hz component of the electromagnetic field and can excite the quasi-BIC mode: a linearly polarized plane wave E=E0eikxy^ (PW) and an azimuthally polarized cylindrical vector beam (AP) [38] (see Supplementary Material). We perform the multipolar decomposition of the incident radiation in terms of vector spherical harmonics (see Supplementary Material) and after numerical integration we obtain the following relations of the magnetic octupolar and magnetic dipolar contributions: |AM(3,0)AM(1,0)|2=0.875 for the PW excitation, and |AM(3,0)AM(1,0)|2=4.854 for AP excitation. Remarkably, an AP vector beam can be used to couple to the BIC mode more efficiently than a linearly polarized plane wave because the AP beam can be decomposed solely to magnetic multipoles with m = 0 that matches the multipolar composition of the modes, while the plane wave consists of electric and magnetic multipoles with different azimuthal indices m. Furthermore, the relative contribution of MO is significantly higher in the AP vector beam than in the plane wave, so that the azimuthally polarized cylindrical vector beam is favorable for excitation of the high-quality resonance in the dielectric nanodisk. Thus, hereafter we focus on the study of SHG in a AlGaAs nanodisk excited by an AP cylindrical vector beam in the vicinity of the BIC point.

2.2 Nonlinear response driven by quasi-BIC state

For noncentrosymmetric materials, the induced nonlinear polarization is defined by the second-order polarizability tensor χ^(2):

where ε0 is the vacuum dielectric constant, E is the electric field inside the particle. For AlGaAs in the principal axis system of the crystal, the tensor of the second-order nonlinear susceptibility contains only off-diagonal elements χijk(2)=χ(2) = 100 pm/V being nonzero if any of two indices i,j,k do not coincide. Thus, the induced nonlinear polarization at the second-harmonic frequency takes the form

We assume that the main axes of the crystalline lattice are oriented along the axis of the Cartesian coordinate system: [100]∥x^,[010]∥y^,[001]∥z^. Thereby, in the case of excitation of the azimuthally polarized modes inside the nanodisk the nonlinear current at the second-harmonic frequency has the form

Here, Eφ(ω) is the electric field distribution inside the particle, which can be approximated by the following expression

where J1 is the first order Bessel function of the first kind, k0=ω/c, ρ=x2+y2 is a radial distance, α1,2 are coefficients of the eigenmodes contributions. Next, using the approximation for the field profile inside the particle, we analyze the multipolar composition of the SH radiation generated by the nonlinear source following the procedure described in Ref. [25]. We deduce that the orbital numbers l of generated multipoles are even for electric and odd for magnetic multipoles, while the azimuthal number takes the values m=±2 only. We find four dominant generated multipoles: aM(2ω)(2,±2), aM(2ω)(4,±2), aE(2ω)(3,±2), aE(2ω)(5,±2). Our analytical calculations show that at k0r∼1 electric multipolar coefficients aE(2ω)(l,±2) rapidly decrease for l≥7, while magnetic amplitudes aM(2ω)(l,±2) decay fast for l≥6. For example, in the vicinity of the BIC point we obtain the following relations in the case of AP beam excitation:

In the case when a nanoantenna supports both a quasi-BIC resonance near the pump frequency and an eigenmode near the SH frequency, we expect a resonant increase in the up-conversion. In addition to the eigenmodes at the fundamental frequency, Figure 4 features the dispersion of the high-quality mode at the second-harmonic frequency lying in the same range of parameters (SH mode). The Q-factor of this mode is large (Q≈200−600) (see Supplementary Material). Both Figures 2 and 4 show two interacting modes responsible for the quasi-BIC formation at the fundamental frequency and selected SH mode which falls in the same parameters’ range and, by symmetry, can be excited by the BIC-driven nonlinear source. Certainly, there are several eigenmodes at the second-harmonic frequency in the considered range of parameters. However, here we specifically focus on the eigenmode that can be addressed by our nonlinear source j(2ω)=2iωP(2ω), namely, only the eigenmode with azimuthal indices m = ±2. This mode is composed of the multipoles that the induced nonlinear current at the quasi-BIC conditions generates. The multipolar decomposition of the SH mode is visualized by circles of different colors. The right inset additionally shows relative contributions of multipoles in the line plot. The SH mode has a more complex multipolar content than the FF modes, and a more sophisticated near-field distribution (top inset), which modifies slightly in the considered parameter range. The radiation pattern of the SH mode is generally multi-lobed and noticeably changes depending on the multipolar composition. The most efficient SHG can be reached at the resonant conditions, close to the dispersion crossing of the quasi-BIC mode and the SH mode.

Figure 4:

Colored circles trace dispersion of the eigenmode of the AlGaAs disk at the second-harmonic frequency (similar to Figure 2, its wavelength is doubled). Sizes of circles scale the total radiative decay (black) and partial multipolar contributions, also plotted as lines in the inset. Black solid lines are dispersion curves of the fundamental-frequency eigenmodes from Figure 2. Insets show the near-field distribution of the electric field magnitude at r/h = 0.729 (top and side views in the upper image) and far-field diagrams at r/h = 0.62, 0.667, 0.729 and 0.775 (the bottom row). E/M – electric/magnetic multipoles, Q – quadrupole, O – octupole, H – hexadecapole, T – triacontadipole.

As a pump, we choose an AP beam in order to maximize the overlaps of the modes with the incident field. We numerically calculate dependencies of multipolar compositions of the linear scattering and SH radiation on the disk radius in the case of AP beam excitation (for comparison with the case of PW excitation see Supplementary Material). The multipolar contents of both FF and SH fields are found to be in excellent agreement with the results of our theory that justifies its validity. Figure 5 shows two-dimensional maps of linear scattering of the AP beam and the SH power. We recover maps of linear scattering, MD and MO contributions in linear scattering on the nanodisk radius by calculating slices at fixed pump wavelengths. We observe strong excitation of a magnetic octupole (MO) in the proximity of the BIC point so that the high-Q mode branch can be clearly traced in the map of linear scattering. Next, we model SHG and obtain a sharp enhancement peak in the SH intensity. The resonant area in the parameter space is visualized with high contrast in Figure 5D. These our results fully explain the recent experimental observations reported in Refs. [39], [40].

Figure 5:

Maps of linear scattering of the AP beam by the AlGaAs nanodisk (A) and the SH generated power (D). Panels (B) and (C) correspond to magnetic dipole and octupole contributions to the linear scattering. Solid pink lines correspond to the dispersion curves of FF eigenmodes, solid blue line is the dispersion curve of SH mode. Pink dot marks a BIC state. Dashed lines are two slices of λ₀ = 1580 nm and λ₀ = 1600 nm. Insets in the panel (D) show electric field distribution at the BIC point (in the pink frame) and at point of intersections of dispersion curves (at maximum SH power) (in the yellow frame).

Figure 6 shows two slices of 2D color maps in Figure 5 for two close wavelengths: the slice at λ₀ = 1580 nm contains the BIC point and the slice at λ₀ = 1600 nm contains the intersection of the SH and FF dispersion curves. If we consider the fixed pump wavelength λ₀ = 1580 nm, we see that MO dominates in the linear scattering and the maxima of the scattering cross-section and the SH power for this slice correspond to the BIC point (r = 457 nm, h = 645 nm). We obtain that in this case magnetic quadrupole (MQ) aM(2ω)(2,±2) dominates the emission, and the SH radiation pattern is mainly quadrupolar. For the wavelength λ₀ = 1600 nm corresponding to the intersection of the FF high-Q mode and SH mode dispersion branches, the SHG enhancement is even higher. In this case major multipoles are aM(2ω)(4,±2), aE(2ω)(5,±2), that matches the multipolar content of the SH mode.

Figure 6:

Multipolar decomposition of linear scattering (A,C) and second-harmonic generation (B,D) spectra for two slices in Figure 5. Top and bottom panels correspond to pump wavelengths λ₀ = 1580 and 1600 nm, respectively. Colored lines show different multipole contributions to the linear scattered and the second-harmonic fields. Insets show far-field diagrams of SH radiation at the BIC point (λ₀ = 1580 nm, r = 457 nm, h = 645 nm) and at the point of dispersion curves intersection (λ₀ = 1600 nm, r = 470 nm, h = 645 nm). E/M – electric/magnetic multipoles, D – dipole, Q – quadrupole, O – octupole, H – hexadecapole, T – triacontadipole.

Because the distribution of the AP beam intensity is inhomogeneous, here we define the scattering efficiency as the ratio of the scattering power at the fundamental frequency P_FF to the energy flux through the geometrical cross section of the particle in the waist plane PAP:

Here, S‾=12Re(E×H∗) is a time-averaged Poynting vector, n∥z^ is a normal vector. The scattering cross section σsca=σ˜S is much larger than the geometric cross section S (S=πr2 for AP and S=2rh for PW) of the disk near the BIC-point especially for the AP beam (σ/πr2≈40) because of the high-quality factor of the MO resonance. In this case, it would be incorrect to define the conversion efficiency as the ratio of the generated power PSH and PAP because these values become comparable near the BIC point. But the SH generated power is much smaller than the scattered power at the fundamental frequency because of the high-quality factor of the excited BIC mode. Thus, to take into account the strong interaction of the incident radiation with the resonant mode, we define the SH conversion efficiency as the ratio of the total SH radiated power PSH to the radiated power at the fundamental frequency PFF:

We obtain that the absolute maximum of SHG corresponds to the intersection point and ideally, it may reach about 1% conversion. Also we calculate the conversion efficiency that doesn’t depend on the incident power

The coupled-mode considerations suggest scaling the SHG efficiency (15) depending on detunings from the resonances at the fundamental and SH frequencies:

Here, QFF and QSH are quality factors of resonances, ωFF and ωSH are resonance frequencies, γFF and γSH are damping constants. The highest yield of SHG can be attained in the double-resonant case.

We compare near-field profiles, far-field diagrams and multipolar compositions of the SH mode and the SH radiation near the peak SHG for two types of incident radiation. We obtain that at the parameters of the maximum SHG, the multipolar composition of the SH generated field, the SH electric field distribution and the SH radiation pattern are almost the same as those of the SH eigenmode, see Figures 4, 5 and 7. This explicitly shows that even though the generated SH radiation is composed of several different multipoles, the only one SH mode characterized by m=±2 is predominantly generated.

Figure 7:

Comparison of the far-field diagrams and multipolar compositions of the SHG electromagnetic field and the disk's eigenmode at the SH frequency for parameters of the maximum SHG efficiency (λ₀=1600 nm, r=470 nm, h=645 nm): 1 - SH mode, 2 - AP beam excitation, 3 - PW excitation. E/M–electric/magnetic multipoles, Q–quadrupole, O–octupole, H–hexadecapole, T–triacontadipole

2.3 Coupled model of SHG with pump depletion

In our numerical modeling, we take into account the nonlinear effects of back-action of the second-harmonic radiation on the fundamental frequency radiation. In the case of highly efficient double-resonant excitation, the electric field amplitudes inside the particle at the fundamental and second-harmonic frequencies may become comparable, and the nonlinear correction of the polarization Pnl(ω) at the fundamental frequency can no longer be neglected. Thus, to simulate the back-action nonlinear effect, we numerically solve the Helmholtz equation for the electric field with nonlinear sources at the fundamental and SH frequencies simultaneously:

The expression for the nonlinear polarization at the fundamental frequency Pnl(ω) in the principle crystalline axis system has the following form:

The nonlinear source at the second-harmonic frequency is given by Eq. (8).

In the conditions of the resonant excitation of the magnetic-octupole quasi-BIC mode and the high-quality mode at the second-harmonic frequency, the electric field can be represented as a superposition

where ϵ1,2(t) are the time-dependent amplitudes, |ϵ1,2|2 are the energies of the modes at the frequencies ω1,2, Δω=ω2−2ω1, |Δω|≪ω1,2 is a detuning from the synchronism. Assuming the amplitudes slowly varying in time, with the characteristic time scale τ=|1ϵ1,2dϵ1,2dt|−1≫2πω1,2, one can write truncated equations

where the coefficient

is proportional to the overlap integral, γ‾1,2 are decay rates of the quasi-BIC and SH modes.

We calculate the scattering and SHG conversion efficiencies taking into account coupling of the excited SH and FF modes and compare them with the numbers obtained in the undepleted pump approximation. The results for two types of incident radiation for parameters corresponding to the maximum SHG efficiency (λ0 = 1600 nm, r = 470 nm, h = 645 nm) are summarized in Table 1.

Values of the pump intensity and the energy flux through the geometric cross section of the disk used in our calculations are given in Table 2. S denotes the geometrical cross section of the particle. In the case of AP beam excitation, Imin and Imax are the minimum and the maximum values of the total field intensity of the AP beam in the focal plane, P is the energy flux through the geometrical cross section of the particle in the focal plane. The characteristics of AP beam are illustrated in Fig. S3 of Supplementary Material. The plane-wave intensity Imin=Imax is spatially constant. Taking into account a finite pulse duration, we have refined the value of the total conversion efficiency ρ˜ for AP beam and compared it with data from experimental works [39], [40] (see Supplementary Material). The calculated total second-harmonic efficiency 9.5×10−5 W⁻¹ for 2 ps lasers pulses is approximately twice higher than the estimate in Ref. [40] based on the experimental data.

We obtain that in the case of AP excitation, nonlinear effects of back-action become noticeable due to efficient excitation of the high-quality modes. Thus, coupling of electromagnetic fields at the fundamental and second-harmonic frequencies should be taken into consideration when the conversion efficiency reaches 1%. In the case of PW excitation, the maximum SHG is also reached for the parameters corresponding to the intersection of the eigenmodes dispersion curves, however, the total generated power turns to be lower than in the case of the AP beam pump. As noted above, it can be explained using multipolar analysis of the incident radiation. The multipolar composition of the incident plane wave is more complicated than multipolar composition of the AP beam, the spectrum of linear scattering contains several multipoles with different values of the azimuthal index m. As a result, the electric field distribution at the BIC point is distorted and appears dependent on the azimuthal angle φ. We still observe a resonance at the BIC point in the linear scattering spectrum but the corresponding scattering efficiency is an order of magnitude lower than in the case of AP excitation. Since the magnetic octupolar relative contribution with m=0 is smaller in the plane wave, the excitation of the quasi-BIC state is less efficient.

In the next section, we consider the BIC-inspired nanoantenna design where the second-harmonic radiation significantly influences the linear scattering. Such nanoresonators can be particularly promising for realization of the efficient frequency downconversion at the nanoscale.

2.4 BIC-inspired nonlinear enhancement in a hybrid metal-dielectric nanoantenna

One type of BICs which has been widely studied is the symmetry-protected BIC occurring at Γ point of the bandstructure of a periodic leaky lattice. In periodic systems, when coupling of a certain resonance to the radiation modes is forbidden by a symmetry mismatch, a symmetry-protected BIC is formed [32], [41]. Such symmetry-protected BICs cannot be excited directly under normal plane wave incidence, while they can be excited externally by breaking the in-plane C2 symmetry of the system, for example, by introducing a defect in the nanodisks to open a radiation channel and transform the ideal BIC to the quasi-BIC with a finite Q-factor [42]. In particular, by properly designing a nanodisk metasurface, it can support the longitudinal magnetic dipole resonance at Γ point at the frequency below the diffraction cut-off. Such MD mode does not couple to the free-space radiation forming a BIC. For isolated nanoresonators and nanoparticle oligomers, the out-of-plane Mie-type magnetic modes can be excited efficiently via an azimuthally polarized beam [43], [44], [45]. For periodic systems, this way is not feasible because the incident radiation in each unit cell has to be a structured field. However, inspired by the symmetry-protected BIC mechanism, here we suggest a novel approach to design and access the high-quality modes in isolated Mie nanoresonators.

A schematic illustration of the proposed principle is shown in Figure 8A. As was described in Ref. [16], when an electric current is placed near a perfect-electric-conductor (PEC) surface, the excited free electron oscillations in the adjacent PEC affect the near-field and far-field properties of such composite system. The effect can be reproduced with an oppositely oriented image of the electric current by using the image dipole model. Using this analogy, a periodic system of MD resonators can be transformed to a single MD resonator surrounded by a PEC-like box. As an example, we employ a dielectric disk surrounded by gold (Au) nanobars, as shown in Figure 8B. Remarkably, by shrinking a periodic system into a single nanoresonator, the quasi-BIC modes can be directly and efficiently excited under structured beam illumination that facilitates nonlinear frequency conversion.

Figure 8:

(A) Top: Schematic of the equivalent transformation between circular polarization currents excited in a periodic lattice of MD resonators and currents excited in a single MD resonator surrounded by PEC structures. Bottom: The transformation from a 2D periodic array of MD resonators to an isolated nanostructure composed of an MD resonator surrounded by gold nanobars. (B) Linear scattering and its multipolar decomposition for the hybrid AlGaAs-disk-with-Gold nanostructure under AP pump illumination. (C) SHG efficiency calculated using the undepleted pump model (solid line) and the coupled nonlinear model with back-action (dashed line). Inset shows the SH radiation pattern at the resonance.

We then take the AlGaAs nanodisk supporting a magnetic dipole resonance near 1550 nm and study the nonlinear performance of our proposed design. The geometric parameters of our structure are determined to be the following: the radius and height of AlGaAs nanodisk are r₀ = 237.5 nm, h0=400 nm; the length, height, and thickness of the gold structure are L=575 nm, h2=2h0, and t=100 nm. We suppose that an azimuthally polarized beam, with the maximum intensity Imax=2.16GW/cm2, focused by an objective with the numerical aperture NA=0.9 is illuminating our sample. Our recent experimental study [16] ensures that the pump intensities ∼3GW/cm2 do not cause damage of gold in the near-infrared range at the 1–2 ps laser pulse duration. Figure 8B shows the numerically computed linear scattering featuring magnetic multipoles with m=0. A narrow dip appears around 1550 nm corresponding to excitation of the high-Q mode. The quality factor of this cavity mode is estimated to be Q=554 which is much larger than that of typical MD resonance supported by an individual nanoresonator with (Q∼ 9) [22]. Its longitudinal nature enables optimal overlapping with an AP pump.

We further examine the SHG process from the designed AlGaAs-Gold hybrid nanoresonator. In order to show the importance of our advanced simulation method introduced in Section 2.3, we compare the calculated results based on the undepleted pump and the coupled back-action models in Figure 8C.

In the undepleted pump model, we follow two steps to simulate the nonlinear response. First, we calculate the linear response at the fundamental wavelength and obtain the nonlinear polarization induced inside the structure based on Eq. (8). We then employ the nonlinear polarization as a source to simulate the electromagnetic response at the harmonic wavelength. In this two-step procedure, we disregard the influence of the harmonic waves on the pump field. As seen in Figure 8C, the SHG emission is significantly enhanced when exciting the designed high-Q resonance at the pump frequency. Near the resonance, we observe a clear difference in the nonlinear emission power when using the coupled back-action model as compared to the undepleted pump model: at the pump wavelength of 1549 nm, the total SHG power is around 0.593 W and 0.546 W for these two cases, respectively. This indicates that neglecting coupling from SH field to the pump field in the undepleted pump model results in more than 8% error in the nonlinear simulation. Under AP pump illumination, the SHG signal can be boosted in our proposed configuration leading to nearly 2000-fold enhancement of the SHG emission power as compared to the case of SHG driven by the conventional MD resonance in a single free-standing nanodisk.