Efficient frequency conversion with geometric phase control in optical metasurfaces

Metasurfaces have appeared as a versatile platform for miniaturized functional nonlinear optics due to their design freedom in tailoring wavefronts. The key factor that limits its application in functional devices is the low conversion efficiency. Recently, dielectric metasurfaces governed by either high-quality factor modes (quasi-bound states in the continuum) or Mie modes, enabling strong light-matter interaction, have become a prolific route to achieve high nonlinear efficiency. Here, we demonstrate both numerically and experimentally an effective way of spatial nonlinear phase control by using the Pancharatnam-Berry phase principle with a high third harmonic conversion efficiency of $10^{-4}$ $1/W^2$. We find that the magnetic Mie resonance appears to be the main contributor to the third harmonic response, while the contribution from the quasi-bound states in the continuum is negligible. This is confirmed by a phenomenological model based on coupled anharmonic oscillators. Besides, our metasurface provides experimentally a high diffraction efficiency (80-90%) in both polarization channels. We show a functional application of our approach by experimentally reconstructing an encoded polarization-multiplexed vortex beam array with different topological charges at the third harmonic frequency with high fidelity. Our approach has the potential viability for future on-chip nonlinear signal processing and wavefront control.


Introduction
Nonlinear metasurfaces provide fine-grained control over the generation and efficient manipulation of new frequencies. Freed from the constraints of classical nonlinear optics [1] , nonlinear metasurfaces feature versatile functionalities by manipulating the amplitude, phase, and polarization of the generated light. Past research led to essential applications like higher harmonic generation and wave mixing [2,3] , beam shaping, wavefront control at higher harmonics [4] , vortex beam generation [5] , multiplexed holography [6] , and nonlinear chirality [7,8] . Many of these applications were demonstrated with plasmonic or hybrid plasmonic-dielectric metasurfaces that are governed by dipole electric or magnetic modes but yield a low conversion efficiency [9,10] in the visible and near-infrared regime. Concerning practical applications, the low conversion efficiency is a severe bottleneck. However, introducing all-dielectric nanostructures made of high refractive index materials with strong nonlinear coefficients allows us to boost the nonlinear conversion efficiency by several orders and introduce certain functionalities [11][12][13][14][15][16][17][18][19][20][21][22] . These high index nanostructures support diverse spatial electromagnetic modes to couple the material nonlinear coefficients to the external electromagnetic field. Two important categories of these modes in the context of nonlinear optics are the Mie modes and the bound states in the continuum (BICs). Mie modes are resonances in dielectric resonators that show a strong electromagnetic response, allowing to confine the local electromagnetic field inside a subwavelength nanostructure and therefore enhance the light-matter interaction. Consequently, Mie resonances can achieve high nonlinear conversion efficiencies. On the other hand, dielectric nanoresonators can also support BICs. Compared to the Mie modes, BICs feature an infinite quality factor because they are modes decoupled from the free space radiation spectrum [23,24] . By perturbing the symmetry of the BICs, it is possible to transform the BICs into a quasi-BICs (QBICs) capable of coupling to the freespace modes [24][25][26] , which further improves the light-matter interaction and leads to high nonlinear conversion efficiencies. [26,27] Note that design of nonlinear metasurfaces often target resonances at the fundamental wavelength to increase the nonlinear conversion efficiency. However, dual resonance techniques are also useful, where one resonance is placed near the fundamental wavelength and another at the generated harmonic wavelength. [28,29] Besides addressing the high nonlinear conversion efficiency, nonlinear wavefront control is another important criterion for functional applications. The most common approach for nonlinear phase-tailoring in metasurfaces is the extension of the Huygens principle to the nonlinear regime.
Here, the phase of the light is manipulated by tuning the shape of a nanoresonator. [13,30] With this approach, to achieve a fine-grained phase resolution, a set of numerous different nanostructures with distinct phase responses need to be designed, which is a complex task and requires heavy computational resources. A more straightforward approach to achieve a fine-grained phase resolution is the geometric phase or Pancharatnam-Berry (PB) principle. Here, the phase of the generated light is adjusted continuously by using the degree of rotation of a single nanostructure.
This approach requires designing and optimizing only a single nanostructure instead of a set of nanostructures for obtaining a high nonlinear conversion efficiency. However, to achieve a high diffraction efficiency and smooth spatial phases, the unit cell size of the metasurface needs to be small. Adversely, a small unit cell size leads to near-field coupling between the nanostructures, which may alter the phase response of the metasurface. Therefore, it is highly challenging to maintain the nonlinear phase relationship resulting from the geometric phase principle, to achieve a continuous phase tailoring over the range of 2π, and at the same time, maintain the desired nonlinear efficiency.
It is to be noted that the previously reported conversion efficiency for third-harmonic generation (THG) in silicon metasurfaces with nonlinear PB phase approach remains very low (~10 −9 1 2 ) compared to the same from the metasurfaces using the Huygens principle (average power conversion efficiency ~10 −4 1 2 ). [12,16,31,32] Often, amorphous silicon is used as a material platform for fabricating metasurfaces to investigate THG due to its mature fabrication technology. The nonlinear PB phase principle requires circularly polarized fundamental excitation compared to the linearly polarized excitation in case of nonlinear Huygens metasurfaces. Therefore, there appears a disadvantage for amorphous silicon in the context of nonlinear PB phase control because of its isotropic nature, as THG from an isotropic material under circularly polarized excitation is forbidden. [33] Consequently, the resonator design with amorphous silicon for nonlinear PB phase control requires special attention for achieving a high nonlinear conversion efficiency for THG. [31,34] In this work, with the goal of achieving a highly efficient nonlinear metasurface together with the PB phase principle, we designed our metasurface based on QBICs that provide high-quality factor modes. Numerically, we achieved a strong third harmonic (TH) signal, mainly originating from the nonlinear enhancement due to the QBIC . The experimental observation of the THG provided   a very high nonlinear conversion efficiency (10 −4 1 2 ) under circularly polarized fundamental excitation. Interestingly, nonlinear measurements indicate that the dominant contribution to the TH response arises from a magnetic Mie resonance as confirmed by an analytical model that is based on the principle of coupled anharmonic oscillators. Further, to establish the functional capability of our nonlinear metasurface, we experimentally reconstructed an encoded vortex beam array with different topological charges at the triple frequency with high fidelity, as illustrated in Figure 1a. We applied a design algorithm optimized to overcome the challenges of optical vortex beam generation and demonstrated the successful encoding of information in multiple channels of the TH signal simultaneously.

Numerical design of nanoresonator and metasurfaces
In our design, we use a two-dimensional array of cylindrical nanoresonators made of amorphous silicon. The unit cell of the metasurface consists of a cylindrical nanoresonator with an off-centered air hole, which alters the rotational symmetry of the cylinder from ∞ to 1 as shown in Figure 1a.
Note that the 1 rotational symmetry helps to introduce the geometric phase to the system, followed by the selection rule for third-harmonic generation under circularly polarized fundamental illumination [34] . The geometrical parameters of the unit cell were optimized based on the QBIC mode to achieve a high nonlinear conversion efficiency for THG under the fundamental excitation of wavelengths around 1300 nm. For this optimization, we use full-vectorial electromagnetic simulations with the finite-element-method implemented in COMSOL Multiphysics (more details can be found in the Supporting Information). The obtained linear transmission of the metasurface for a period of 664 nm under circularly polarized light is shown in Figure 1b. One can observe a broad resonance dip at a wavelength of 1300 nm with a linewidth of ~70 nm, resembling a magnetic Mie mode, and a sharp peak at a wavelength of 1325 nm with a linewidth of ~8 nm, which corresponds to QBICs. The magnetic field plot at the wavelength of 1300 nm (top inset of Figure 1b) shows the magnetic dipole nature of the resonance located in the x-y plane. Figure 1c illustrates the corresponding electric field plots in the x-z and y-z planes for the same magnetic field in the x-y plane (repeated) at 1300 nm. The electric field forms a vortex in the x-z and y-z planes, which is typical for a magnetic-dipole Mie resonance with in-plane magnetic dipoles, in particular aligned with the x-axis and y-axis. [35] The magnetic dipole mode has the ability to enhance the nonlinear conversion efficiency as they allow for efficient coupling of an electromagnetic field to the nonlinear susceptibility of a material. In contrast, the embedded sharp peak at 1325 nm shows an entirely different mode profile compared to the broad resonance at 1300 nm. At 1325 nm, the electric field plot shows an asymmetric vortex in the x-y plane (bottom inset of Figure 1b), which is enhanced in the vicinity of the air hole. The corresponding magnetic field in the y-z plane is parallel to the cylinder axis, as shown in the lower row of the Figure 1d, along with the same electric field in the x-y plane.
To better understand the QBIC, we consider the nanoresonator without the small air hole. The electric and magnetic field plots of such a cylinder are shown in the upper row of Figure 1d. As visible in the magnetic field plot, this geometry supports a vertical magnetic dipole mode, parallel to the cylinder axis; therefore, the electric field of the mode forms a symmetrical vortex in the x-y plane. The mode is incapable of coupling to the free space as the mode's symmetry in the x-y plane is mismatched to the symmetry of a plane wave at normal incidence, and the unit cell subwavelength period prevents any wavevector which is not normal to the metasurface. [36] These properties eliminate possible channels for radiative decay, resulting in a high-quality factor and a long lifetime. The long lifetime of these so-called BICs is vital for high nonlinear conversion efficiency as it allows for a prolonged interaction between the trapped energy and the nonlinear susceptibility of the material. [25] By incorporating an off-center hole within the cylinder, one can break the in-plane symmetry of the nanoresonator, open radiation channels, and the BIC is then transformed into quasi-bound states in the continuum (QBICs), which is capable of free-space coupling. The corresponding electric and magnetic fields for the QBIC for an array of the cylinders with holes of radius 40 nm are shown in the lower row of the Figure 1d. Compared to the electric field plots for a cylinder without a hole, the hole leaves the electric field vortex mostly intact but results in a field enhancement in its vicinity. This property leads to a nonzero electric dipole moment in the x-y plane, compatible with the symmetry of the plane wave at normal incidence.
As mentioned above, one can use the hole in the cylinder to introduce the Pancharatnam-Berry (PB) phase for circularly polarized light. The qualitative relationship between the rotation of a nanostructure with 1 rotational symmetry with respect to the laboratory frame and the obtained phase for THG is given by [31,37] : where, and are the acquired nonlinear phase shifts when the fundamental and TH light has the same (co-), and the opposite (cross-) circular polarization, respectively, and = ±1 defines the handedness of the circular polarization state. The principle is illustrated in Figure 2a. For functional applications such as in nonlinear holography, a high spatial density of nanoresonators is necessary to achieve high-quality holographic images as the image quality increases with the spatial phase resolution. On the other hand, increasing the density of resonators increases the near-field coupling between the adjacent unit cells [38] , leading to a deviation in the TH phase relationship as presented in the Equations 1 and 2. It results in an increased crosstalk between the different polarization channels and introduces a phase noise, limiting the diffraction efficiency. Therefore, due to near-field coupling, there always exists a trade-off between the PB-phase relationship and a high spatial phase resolution To assess the trade-off between the near-field coupling and the phase, we numerically investigated the TH intensity and the corresponding nonlinear phase in our design for different periods. Figure   2b shows the color plot of the transmission spectra within the wavelength range of 1100-1500 nm as a function of the unit cell period, (500-900 nm). The transmission spectrum presents a sharp peak within a much broader resonance dip, corresponding to the QBIC as indicated previously in angle of a single nanoresonator between 0° and 180° changes the phase of the generated TH light within 0 to 2 for co-polarization and within 0 to 4 for the cross-polarization due to the different dependency of the acquired nonlinear phase on the rotation angle of the nanoresonator (2 and 4 in co-and cross-polarizations, respectively). Besides, the nonlinear phase relationship follows higher linearity with larger periods, as the near field coupling becomes weaker with the larger period. Therefore, the unit cell size is crucial for optimizing the system for nonlinear PB phase control, while maintaining a uniform TH intensity.

Experimental outcomes and analytical modelling
Based on our design, we fabricated 10 different metasurfaces of size 100x100 µm² for five different   Performing a Circle Hough Transformation on several SEM images, we obtained the average radius of the cylinders to be (182.5 ± 3.5) nm, which matches to our design (182.5 nm). However, the hole radius after the fabrication appears to be smaller than the design and increases with the period due to the decrease of the proximity effect during electron beam lithography. Therefore, instead of having a uniform radius of 40 nm, the hole radius in our sample varies from 20 nm for a period of 562 nm to 35 nm for 811 nm. The resonance dip for the horizontally polarized excitation lies at a shorter wavelength than the same for the vertical polarization for the same period. In agreement with the simulation, the resonance frequency shifts to longer wavelengths with an increasing period. Moreover, the resonance shape differs for both polarizations; the dips for the horizontal polarization are symmetric, while it appears as asymmetric for vertical polarization, which may result from the influence of the QBIC mode. However, the anticipated peak in the transmission spectra due to QBIC as observed in the simulation (Figure 1d) is not observed experimentally.  Also note that the estimated efficiency is based only on the zeroth diffraction order of the twodimensional metagrating. As the TH wavelength is smaller than the period of the unit cell one can expects higher diffraction orders to appear. The higher diffraction orders appear with an angle greater than 45° and thus are not collected by the microscope objective, as the numerical aperture only allows for collecting light within a cone of 24°. Therefore, if we include the 9 diffraction orders (because of the two-dimensional grating), the overall TH conversion efficiency will be one order of magnitude higher. The efficiency remains the same in the case of the corresponding phase gradient metasurface. Compared to an unstructured silicon film of the same thickness, the patterned metasurface shows a 10,000 times higher TH intensity (see Supporting Information). It is worthy to mention that the THG peak power conversion efficiency (in the zeroth order) of the metasurface, when illuminated by a vertically polarized fundamental beam is T HG,lin =̂T HG (̂i n ) 3 ≈ 10 −13 1 2 (10 −3 1 2 is the average power conversion efficiency). By using a circularly polarized fundamental beam, we lose 50% of the power during excitation, leading to ~( in the vicinity of that peak, our fabricated metasurface presents a broadband nonlinear response in the experiment. For a better assessment of this discrepancy between the experimental and the numerical outcomes, we adopted a model based on a coupled anharmonic oscillator [39] . In this model, we define the amplitude of the nanostructure's magnetic Mie mode as which can interact with an external excitation and couple to the free space. We also define the BIC amplitude as which cannot be coupled directly to the free space but is coupled to . Therefore, one can obtain two coupled differential equations as follows: Where / ( ) / , / represent the amplitude, damping, and resonance frequency of the Mie mode (D) and the BICs (B), respectively, and is the coupling constant between the oscillators.
Further, only the Mie resonance is driven by an external electromagnetic field 1 , since the BIC is a dark mode. The anharmonic term ]. [39] The nonlinear response of the metasurfaces can be estimated by using the extracted values of / , / and from the fitting of the extinction spectra experimentally obtained from the transmission spectra shown in Figure 3c. The fit values can be found in Table 1  properties are mostly remain unchanged. Additionally, the nonlinear susceptibility of the material is also not considered in this model but only the nonlinear behavior based on the resonance itself, which also influence the deviation between the model and the experimental data. Besides, the coupled anharmonic oscillator model, temporal coupled mode theory (CMT) could also be useful to describe and analyze the third harmonic response from resonant nonlinear systems [40] . However, our model reasonably predicts the overall nonlinear response of the metasurface, and we can separate the individual contributions from the two individual modes. Figure   For further investigation, we fit our model to the numerical transmission spectra as displayed in Figure 1b. The fit values can be found in Table 1. The resonance wavelengths corresponding to the frequencies and of the magnetic Mie mode and the QBIC fall at the minimum of the broad dip and the sharp peak in the transmission spectra, respectively. In contrast to the linewidths for BIC and Mie modes observed in the experiment, the linewidths extracted from the simulation are appeared to be smaller. The fitting presented in Table 1 reveals a linewidth of ~40 nm for the Mie mode, which is close to the value obtained from the experimentally measured transmission spectra. However, in Table 1 the linewidth of the mode is estimated as 0.1 nm, which corresponds to a quality factor of 13.200, agrees with the assumption that the dark mode in this resonator system resembles a BIC mode. In addition, the coupling constant can be estimated as = 16 nm, which is smaller than the linewidth of the magnetic Mie mode, indicates a connection with the Fano resonant regime associated with weak coupling [41] . Therefore, in the numerical case, our design supports BIC and the coupled anharmonic oscillator established by equation 3 and 4 describes a QBIC for the set of parameters presented in Table 1. Figure Table 1. Therefore, the linewidth of the QBIC is the key parameter that makes our experimental observations significantly different from the numerical results, the linewidth of the QBIC appears to be around 12 nm in the experiment. Further, the model parameters can be extracted from the simulated transmission spectrum shown in Figure 1d Figure 1d becomes less pronounced and resembles the measured transmission in Figure   3c, as the linewidth is increased from 0.1 nm to 6.5 nm. Besides the linewidth, the spectral distance of the BICs and Mie resonance appears to influence the nonlinear response of the system as well.
In total, the dominant contribution in the experiment stems from the magnetic Mie mode and not from the QBIC since the linewidth of the QBIC, which represents the loss of this mode, appeared to be too large to drive the QBIC effectively. Reason for the large losses, which limit the contribution of the QBIC mode can be mostly fabrication related. This includes geometrical variations between the resonators during the lithographic patterning, as explained in Ref. [42] , and changes in the refractive index during the reactive ion etching and the other fabrication steps. In addition, potential surface roughness of the resonators, which also leads to a broadening, was not included in our design. Overall, these parameters can all limit the quality factor of the mode. Note that QBIC in silicon metasurfaces with high quality factors were demonstrated. [43,44]  To realize the nonlinear phase-tailoring properties experimentally, we investigated THG from our fabricated metasurfaces with the desired phase gradient based on the PB phase principle for different unit cell periods, leading to a beam deflection (for details, see Supporting Information).
Depending on the period, one expects diffraction orders at around 4.64° (co-polarization) and 9.29° (cross-polarization) for a TH wavelength of 430 nm. As the phase gradient resembles a blazed grating, the diffraction spots can be assigned to different phase gradients resulting from the phase factors of the PB-phase (±2 ↔ ±4 , ±1 ↔ ±2 ), as explained in the Supporting Information. We spatially resolve the TH signal in the Fourier-space by imaging the back focal plane of the microscope objective that collects the light from the metasurface. The measured TH intensity in different diffraction orders, along with the estimated diffraction efficiencies for coand cross-polarization, are shown in Figure 5. The k-space image for the co-polarized TH light shows a prominent spot at 1 st diffraction order (at the desired angle) with unwanted weaker contributions at the 0 th and 2 nd order location. In contrast, the cross-polarized TH light appears mainly in the 2 nd diffraction order. The 0 th order diffraction spot in the center of the Fourier space does not carry any PB phase. We calculated the diffraction efficiency for a particular diffraction order as the quotient of the intensity of that diffraction order (marked by the dotted white lines) to the overall TH intensity.
In co-polarization, the first diffraction order (proportional to 2 ) appears with very high diffraction efficiency (~90%) for periods ranging from 664 nm to 811 nm. The high diffraction efficiency is maintained for a broad range of fundamental wavelengths. However, for smaller periods (609 nm and 562 nm), the diffraction efficiency decreases dramatically due to the upsurge of the near-field coupling strength between the adjacent nanoresonators. Moreover, the diffraction efficiency for the −1 (Figure 5b) and −2 (Figure 5a) diffraction orders are negligible as the phase gradient resembles blazed gratings. The same trend of variation of the diffraction efficiency is also observed for cross-polarized detection configuration. For cross-polarization, the phase shift is proportional to 4 and therefore, the maximum diffraction efficiency occurs at the +2 diffraction order ( Figure 5e). Here, we achieve a maximum diffraction efficiency of ~80%. However, the overall diffraction efficiency is estimated as lower for the cross-polarization compared to the co-polarized configuration. This is primarily because the TH phase in the cross-polarization changes with 4 while in the co-polarization case the phase changes with 2 , decreasing the phase resolution by a factor of 2. Consequently, the diffraction efficiency decreases for the cross-polarization compared to the co-polarization.
To establish the functional capacity of our metasurface, we designed a vortex beam array for its realization at the triple frequency. The design was based on the spatial distribution of the TH phase to accomplish the desired nonlinear diffractive meta-element. We chose a 3x3 arrangement of This effect is stronger for the vortex beam with low topological charges, e.g., for m = -2, as the intensity is concentrated over the smaller area on the camera sensor.

Conclusion
In this work, we demonstrate both numerically and experimentally a nonlinear silicon metasurface design. [13,16] In contrast, the nonlinear extension of the PB phase principle, which requires a more straightforward design and a better-simplified nanofabrication process, can end up with very high conversion and diffraction efficiencies even at a circularly polarized fundamental excitation, along with a continuous and full control of the nonlinear phase, as evidenced by this work. Finally, we established the functional capability of our metasurface by illustrating the experimental reconstruction of a polarization multiplexed encoded vortex beam array supported by varying topological charges, with high fidelity and high conversion efficiency.
The concept provides a pathway to realize the miniaturized control of the high conversion and diffraction efficiencies to be applicable in nonlinear nanophotonic devices where advanced phase control with low spatial footprints is necessary, such as nonlinear optical information processing and wavefront engineering.

Optical experiments
The linear transmission spectra were measured using a white light Laser source (Fianium Whitelase) in the spectral range of 1150 nm to 1300 nm. Further, we controlled the input polarization by a linear polarizer. A 10x infinity-corrected microscope objective collected the metasurface's transmission. The light was routed to a spectrometer (Andor Shamrock Kymera 193) with a set of lenses. The spectrometer was equipped with an InGaAs detector (Andor IDus 491A InGaAs) suitable for measurements in the infrared. The spectral resolution of the transmission setup is 2.43 nm.
We use an Optical Parametric Oscillator with a pulse length of 200 fs and an 80 MHz repetition rate as a laser source for the nonlinear optical experiments. We change the wavelength between 1200 nm and 1350 nm during the experiment but keeping the average laser power at 50 mW. The laser light is focused on the metasurface by a 50 mm achromatic lens optimized for the infrared resulting in a beam waist size of approximately 45 µm. Therefore, the average pump peak power density of 0.32 GW cm 2 . To collect the THG light from the metasurface, we use a 50x microscope objective with a numerical aperture of 0.42. By a set of two lenses, we capture the back focal plane into an imaging spectrometer. The imaging spectrometer (Andor Shamrock 303 equipped with an iDUS 420 detector) allows us to resolve the frequency-converted light's wavelength and different diffraction orders orthogonal to the diffraction grating. Using a combination of linear polarizers and quarterwave plates, we prepare and analyze the circular polarization state. We have used short pass filters to suppress the fundamental Laser beam, to minimize nonlinear contributions from other optical components. The holographic images were captured with a similar setup as discussed above, but an sCMOS camera (Andor Zyla 4.2) replaces the spectrometer to capture a high-fidelity image of the whole Fourier space.

Vortex beam array design
The vortex beam array metasurface design exploits the same QBIC metasurface as discussed above, but a computer-generated holographic (CGH) algorithm determined every nanoresonator's rotation on the metasurface. This algorithm is optimized for vortex beam generation at the Fourier plane. There is a particular need for optimized algorithms since the generation of optical vortex beam arrays has several challenges, such as low quality and an unfavorable energy distribution among the vortex beams. Further, the algorithm allows for polarization and wavelength multiplexing, satisfying the demand for compact, highly functional meta-devices.
To establish the functional capacity of our nonlinear metasurface, we designed a nonlinear Dammann vortex beam arrays for its realization at two circular polarization channels of THG and third channel in linear case. Since for such C1 symmetric meta-atom, one can predict 2 ( = −1) for opposite handedness at the fundamental frequency, and 2 ( = 1), 4 ( = −1) for same and opposite handedness at the THG frequency, respectively. The transmission function can be expressed by Fourier series as follows: Where the inset topological charges are set as 0 = 2, = 2, = 6, n represents the harmonic generation order, which can be 1 (for fundamental frequency) and 3 (for THG). While (n-1) and (n+1) correspond to the same handedness and opposite handedness of the output light with respect to the incident one, respectively. are the optimizing target to get desired diffraction orders as 3x3. As demonstrated, the topological charges of each diffraction order obey the rule of 0 + + for 2 case and 2 0 + + for 4 case. For each diffraction order, the numerical aperture NA = √ 2 + 2 , which changes with the generated frequencies.
The first item in the series expansion predict the location of each diffraction order, which also relates to the generated frequency. Afterwards, by using genetic algorithm method, the optimized rotation angles of the nanoholes within the cylinders can simultaneously obtain the ideal target orders in both linear and nonlinear channels.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.