Multipolar second-harmonic generation from high-$Q$ quasi-BIC states in nanoresonators

Multimode interference and multipolar interplay govern functionalities of optical nanoresonators and nonlinear nanoantennas. Recently, excitation of the high-quality supercavity modes (quasi-BIC states) in individual subwavelength dielectric particles has been predicted to boost the nonlinear frequency conversion at the nanoscale. Here, we put forward the multipolar model which captures the physics behind linear and nonlinear response driven by such high-$Q$ modes in nanoresonators. We show that formation of the quasi-BIC state in the AlGaAs nanodisk can be understood through multipolar transformations of coupled leaky modes. In particular, the hybridized axially symmetric TE-polarized modes can be viewed as superpositions of multipoles, with a basis of parent multipoles constituted mainly by magnetic dipoles and octupole. The quasi-BIC point in the parameter space appears where dipolar losses are totally suppressed. The efficient optical coupling to this state is implemented via azimuthally polarized beam illumination matching its multipolar origin. We establish a one-to-one correspondence between the standard non Hermitian coupled-mode theory and multipolar models that enables transparent interpretation of scattering characteristics. Using our approach, we derive the multipolar composition of the generated second-harmonic radiation from the AlGaAs nanodisk and validate it with full-wave numerical simulations. Back-action of the second-harmonic radiation onto the fundamental frequency is taken into account in the coupled nonlinear model with pump depletion. A hybrid metal-dielectric nanoantenna is proposed to augment the conversion efficiency up to tens of per cent and actualize the nonlinear parametric downconversion. Our findings delineate the in-depth conceptual framework and novel promising strategies in the design of functional elements for nonlinear nanophotonics applications.

Methods typically employed to describe nonlinear harmonic generation at the nanoscale are based on mul-tipolar decomposition of the fields in spherical multipoles [4,[24][25][26][27]. This approach provides a transparent interpretation for the measurable far-field characteristics, such as the conversion efficiency and radiation patterns, since they are essentially determined by the interference of dominating multipolar modes [4]. The nonlinear response depends on both incident polarisation and symmetry of the a specific material or composite structures, such as nanoparticle oligomers [27][28][29]. For instance, the nonlinear response of disk-shaped nanoantennas made of noncentrosymmetric III-V semiconductors (GaAs or AlGaAs) grown along (100), (110) and (111) crystallographic directions exhibits four-fold, two-fold and continuous rotational symmetries, respectively [26,30]. Specifics of the radiation characteristics in each case can be explained by different parities of nonlinearly generated multipoles, so-called nonlinear multipolar interference. This represents an example where multipolar analysis has been proven to be a useful instrument for design of directional nonlinear light emission [31].
A recently suggested approach to trap light in individual subwavelength dielectric nanoresonators is based on high-quality supercavity modes associated with the physics of quasi bound states in the continuum (quasi-BIC) [32][33][34]. These states are in some sense similar to BICs in infinite periodic dielectric structures: their high finesse is due to the destructive interference of several far-field radiation channels. Two interpretations of quasi-BIC formation in the open resonator, based on (i) leaky modes supported by the particle [33] and (ii) multipoles of the electromagnetic field [35,36], were discussed, however no direct correspondence between these two models has been established up to now. Lately, it has been proposed to utilize these high-quality states to enhance the classical nonlinear process of second-harmonic generation in dielectric nanodisks [37].
In this paper, we develop a comprehensive multipolar theory of the second-harmonic generation (SHG) from high-quality quasi-BIC states in a AlGaAs cylindrical nanoresonator (nanodisk). We show that the formation of quasi-BIC states can be naturally understood through multipolar transformations of coupled leaky axially symmetric modes supported by the nanodisk. The strong SHG can be expected in the case of efficient excitation of the mode with the high quality factor. To achieve efficient coupling to the mode at the pump frequency, multipolar composition of the pump source should match the multipolar structure of this mode. We show that azimuthally polarized beam can be used to couple to the quasi-BIC state most efficiently because it contains magnetic multipoles with zero azimuthal numbers only that matches the quasi-BIC multipolar origin and maximizes the modal overlap. Analysing spatial distribution of the induced nonlinear source with the approach put forward in our earlier works [25,26], we reveal that complex multipolar composition of the second-harmonic radiation coincides with the multipolar composition of a particular high-quality eigenmode of the disk at double frequency. We obtain that the efficiency of SHG is strongly enhanced up to several per cent, provided the generated frequency matches the supported resonance, for parameters of the nanoresonator corresponding to the nearly resonant excitation of the quasi-BIC state at the fundamental frequency. We perform full-wave numerical simulations of SHG taking into account nonlinear effects of back-action and propose the BIC-inspired design of a hybrid metaldielectric nanoantenna where the effect of pump depletion is further increased suggesting a promising application for the frequency downconversion.

II. MULTIPOLAR MODEL OF QUASI-BIC FORMATION
We consider a high-index cylindrical dielectric resonator which 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 homogeneous host medium with ε h = 1. The harmonic time dependence of the fields in the form e iωt is implied.
Here we focus on rotationally symmetric TE-polarized modes in the cylindrical coordinate system (E ϕ , H ρ , H z ), FIG. 1. Schematic of the geometry. Azimuthally polarized cylindrical vector beam of frequency ω excites an axially symmetric supercavity mode in the nanodisk. As a result of nonlinear interaction, the second-harmonic light of 2ω is generated. whose electromagnetic field does not depend on azimuthal angle ϕ. With the use of 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/λ 0 . The results of calculations for the AlGaAs nanodisk with fixed height h = 645 nm are summarized in Fig. 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 in the following.
For a nonspherical shape of the resonator the interacting modes can be in general viewed as superposition of the spherical multipoles distinguished by orbital l and azimuthal m indices. Multipolar analysis of the considered TE-polarized azimuthally symmetric eigenmodes suggests that a basis of parent multipoles is constituted mainly by longitudinal magnetic dipoles (MD) l = 1, m = 0 and magnetic octupole (MO) l = 3, m = 0 so that each mode has two multipolar radiation channels, dipolar and octupolar. In Fig. 2, sizes of circles illustrate their relative contributions in multipolar expansions of the modes. The occurrence of high-Q supercavity mode is accompanied by increasing the order of a dominating multipole from l = 1 (MD) to l = 3 (MO) and corresponds to the decoupled magnetic octupole. At quasi-BIC condition (r/h = 0.71 and r/λ 0 = 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 rotationally symmetric magnetic octupole with three lobes. In addition, the black line in Fig. 2 shows dispersion of the mode near the doubled 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 [38].
We then corroborate our rigorous numerical results by an analytical model setting the correspondence of the non-Hermitian coupled-mode theory and multipolar structure. Coupling of two leaky eigenmodes shown in Fig. 2 can be described by the Hamiltonian [33,34] being a square matrix, where E 1 and E 2 are complex frequencies of the modes 1 and 2, W and V are complex coupling parameters including interaction between the modes via free space, W = V * . A scheme of this phenomenological two-level model is shown in Fig. 3(a).  The eigenfrequencies depend on the aspect ratio ξ ≡ r/h and can be deviated from the BIC-point at ξ * = 0.71: . Also the frequency of incident radiation ω s may differ from the resonance frequency ω s =ω + Ω. For our system, we derive the dynamic equations for two modes with slowly varying amplitudes B 1 and B 2 excited by the external pump source in the following form: where the matrixĤ 2 explicitly readŝ Here d 1,2 are the effective dipole moments, and d 01,02 are effective octupole moments in the modes, S d and S 0 are dipolar and octupolar contributions in the incident radiation respectively. The d coefficients govern both the radiative decay of the states and their coupling to each other and to the external field. Using the scattering-matrix method [35], coupledmode equations can be mapped onto the three-state model schematically illustrated in Fig. 3(b). The basis includes two magnetic-dipole states |1 , |2 with m = 0 and the magnetic-octupole state |0 with m = 0. The evolution equations recast as with the amplitude column-vector a = [a 1 , a 2 , a 3 ], and the external source The non-Hermitian Hamiltonian of the structure in Eq. (4) assumes the form where the condition |ω 0 (ξ * )| > d 2 1,2 > d 2 0 is satisfied, and parameters are expressed as follows The τ coefficients describe the hybridization of the dipole and octupole modes. Ohmic losses are neglected.
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 octupolar S 0 and dipolar S d components and their fractions determine the resultant response. We consider two types of pump radiation which carry H z component of the electromagnetic field and can excite the quasi-BIC mode: a linearly polarized plane wave E = E 0 e ikxŷ (PW) and azimuthally polarized cylindrical vector beam (AP) [39] (see Supporting Information). We perform multipolar decomposition of the incident radiation in terms of vector spherical harmonics (see Supporting Information) and after numerical integration we obtain the following relations of the magnetic octupolar Remarkably, AP vector beam can be used to couple to the BIC mode more efficiently than linearly polarized plane wave because 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 AP cylindrical vector beam in the vicinity of the BIC-point.

III. 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 secondorder nonlinear susceptibility contains only off-diagonal elements χ (2) ijk = χ (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 crystalline lattice are oriented along the Cartesian coordinate system: [100] x, [010] ŷ, [001] ẑ. Thereby, in the case of excitation of azimuthally polarized modes inside the nanodisk the nonlinear current at double frequency has the form Here E (ω) ϕ is the electric field distribution inside the particle, which can be approximated by the following expression where J 1 is the first order Bessel function of the first kind, k 0 = ω/c, ρ = x 2 + y 2 is a radial distance, α 1,2 are coefficients of the eigenmodes contributions. Next, using approximation for the field profile inside the particle we analyse 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: a  ∼ 0.08. (11) In the case when nanoantenna supports a quasi-BIC resonance near the pump frequency and another eigenmode near the SH frequency we expect a resonant increase in the up-conversion. In addition to two eigenmodes at the fundamental frequency, Fig. 4 features the dispersion of the high-quality mode at double frequency lying in the same parameters range (SH mode). The Qfactor of this mode is large (Q ≈ 200 − 600) (see SI). Certainly, there are several eigenmodes at double frequency in the considered parameter ranges. However, here we specifically focus on the eigenmode that can be excited by our nonlinear source j (2ω) = 2iωP, 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 SH mode is visualized by circles of different colors. The inset additionally shows relative contributions of multipoles in the line plot. SH mode has a more complex multipolar content than FF modes, and more sophisticated near-field distribution, which modifies slightly in the considered parameter range. The radiation pattern is generally multi-lobed and noticeably changes depending on the multipolar composition. The most efficient SH generation can be reached at the resonant conditions, close to crossing of BIC mode and SH mode.
As a pump, we choose AP beam in order to maximize overlaps of the modes with the incident field. We numerically calculate dependencies of multipolar composition of linear scattering and SH generated field on the disk radius in the case of AP beam excitation (for comparison with the case of PW excitation see Supplementary  Information). The multipolar contents of both FF and SH fields appear essentially the same as we expected from the theory that justifies its validity. Figure 5 shows twodimensional maps of linear scattering of AP beam and 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 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 SH intensity. The resonant area in the parameter space is visualized with a high contrast in Fig. 5(d). These our results fully explain the recent experimental observations reported in Refs. [40,41].  M (2, ±2) dominates in the emission, and the SH radiation pattern is mainly quadrupolar. For the wavelength λ 0 =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 a Because the distribution of the AP beam intensity is inhomogeneous, here we define the scattering efficiency as the ratio of the scattering power at fundamental frequency P FF to the energy flux through the geometrical cross section of the particle in the waist plane P AP : Here, S = 1 2 Re (E × H * ) is a time-averaged Poynting vector, n ẑ is a normal vector. The scattering cross section σ sca =σS is much larger than the geometric cross section S (S = πr 2 for AP and S = 2rh for PW) of the disk near the BIC-point especially for AP beam (σ/πr 2 ≈ 40) because of the high quality-factor of MO resonance. In this case it would be incorrect to define conversion efficiency as the ratio of generated power P SH and P AP because these values become comparable near the BIC point. But SH generated power is much smaller than scattered power at fundamental frequency because of large quality factor of excited BIC mode. Thus, to take into account the strong interaction of the incident radiation with the resonance mode we define the secondharmonic conversion efficiency as the ratio of the total SH radiated power P SH to the radiated power at the fundamental frequency P FF : We obtain that the absolute maximum of SHG corresponds to intersection point and ideally, it may reach about 1 % conversion. Also we calculate conversion efficiency that doesn't depend on the incident power The simple coupled-mode considerations suggest scaling the SHG efficiency (15) depending on detunings from the resonances at the fundamental and SH frequencies: (16) Here, Q FF and Q SH 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 SH mode and SH radiation near the peak SHG for two types of the incident radiation. We obtain that at the parameters of the maximum SHG, the multipolar composition of the SH generated field, electric field distribution and radiation pattern are almost the same as those of the SH eigenmode, see Figs. 4, 5, 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.

IV. 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 doubled frequencies may become comparable, and the nonlinear correction of the polarization P (ω) nl at the fundamental frequency can no longer be neglected. Thus, to simulate the back-action nonlinear effect, we numerically solve the Helmholtz equation for electric fields with nonlinear sources at the fundamental and SH fre-quencies simultaneously: The expression for nonlinear polarisation at fundamental frequency P (ω) nl in the principle crystalline axis system has the following form: The nonlinear source at the second-harmonic frequency is given by Eq. (8).
At the conditions of the resonant excitation of the magnetic-octupole quasi-BIC mode and 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 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 scattering and SHG conversion efficiencies taking into account coupling of 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 I  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 fundamental and second harmonic frequencies should be taken into consideration when the conversion efficiency reaches 1%. For the 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 AP beam pump. As we noticed 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 AP beam, the spectrum of linear scattering contains several multipoles with different values of azimuthal index m. As a result, the electric field distribution at the BIC point is distorted and appears dependent on azimuthal angle ϕ. We still observe a resonance at the BIC point in the linear scattering spectrum but the quality factor of this resonance is lower, and the corresponding scattering efficiency is an order of magnitude lower than in the case of AP excitation. Also, as we have shown, the magnetic octupolar relative contribution with m = 0 is smaller in the plane wave, thereby 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.

V. BIC-INSPIRED NONLINEAR ENHANCEMENT THROUGH HIGH-QUALITY MD RESONANCE VIA AP BEAM ILLUMINATION
One type of BICs which has been widely studied is the symmetry-protected BICs existing at Γ point of a periodic system. In periodic system, when the coupling of a certain resonance to the radiation modes are forbidden by symmetry mismatch, a symmetry-protected BIC is formed [32]. Such symmetry-protected BICs cannot be excited directly under normal plane wave incidence, while they can be excited externally by breaking the in-plane C 2 symmetry of the system, for example, by introducing a defect in the nanodisks to open a radiation channel and transform the ideal-BIC to quasi-BIC with finite Q factor. In particular, by properly designing a nanodisk metasurface, it can support resonant longitudinal mag- netic dipole resonance at Γ point at the frequency below the diffraction limit. Such MD mode does not couple to the free-space radiation channel due to the symmetry mismatch, forming a BIC. For single nanoresonators, the out-of-plane Mie-type magnetic modes can be excited efficiently via azimuthally-polarized beam [42], whereas for periodic systems this is not possible because the input field for each unit cell has to be a structured field distribution. However, inspired by this symmetry-protected BIC mechanism, here we suggest a novel approach to design high-quality mode in isolated Mie nanoresonators by lifting the restriction on the periodic boundary conditions.
A schematic illustration of the proposed principle is shown in Fig. 8(a). As was described in Ref. [16], when an electric or polarization current is placed near a 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 as a single MD resonator surrounded by a PEC-like box.
As an example, we employ a dielectric disk surrounded by Au nanobars, as shown in Fig. 8(b). Remarkably, by shrinking the periodic system into a single nanoresonator, the quasi-BIC MD mode can be directly and efficiently excited under structured beam illumination that further facilitates nonlinear frequency conversion.
We then take the AlGaAs nanodisk supporting a magnetic dipole resonance near 1550 nm as 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 Al-GaAs nanodisk are r 0 =237.5 nm, h 0 =400 nm; the length, height and thickness of the gold structure are L =575 nm, h 2 = 2h 0 and t =100 nm. We suppose an azimuthally-polarized beam, with the maximum intensity I max = 2.16 GW/cm 2 , focused by an objective with the numerical aperture NA= 0.9 is illuminating our sample. Figure 8(b) shows the calculated linear scattering.
A narrow dip appears around 1550 nm corresponding to excitation of the high-Q quasi-symmetry-protected-BIC MD mode. The quality-factor of this supercavity mode is estimated to be Q = 554 which is much larger than the typical MD resonance supported by an individual nanoresonator (∼ 9) [22]. Its longitudinal MD nature enables optimal overlapping with a 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 Sec. IV, we compare the calculated results based on the undepleted pump model and the coupled back-action model.
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 the formula (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. Figures 8(c) and (d) show the calculated SHG power and the corresponding nonlinear multipolar contributions for different pump wavelengths using the undepleted pump model and the coupled nonlinear model with back-action, respectively. As can be seen, the SHG emission can be significantly enhanced when exciting the designed high-Q resonance at the pump frequency. 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 significantly in our proposed configuration leading to nearly 2,000-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.

VI. CONCLUSION
We have explicated the multipolar nature of quasi-BIC states in subwavelength dielectric resonators supporting high-quality supercavity modes called quasi-bound states in the continuum and discussed their applications for nonlinear nanophotonics. Using our multipolar model, we have analyzed optimal conditions for the efficient excitation of quasi-BIC states in high-index dielectric nanodisks under structured light illumination. In particular, we have explained the multifold increase of the secondharmonic conversion efficiency for the case of azimuthally polarized cylindrical vector beam illumination compared to the linearly polarized plane wave excitation. Implementing numerically the coupled nonlinear model, we have clarified values of SHG efficiencies and discussed pump depletion in isolated high-quality nanoresonators.