Chirality is a fundamental physical property in nature that describes that an object cannot completely coincide with its mirror image through any continuous rotation and translation operation [1–3]. This unique asymmetry arises universally across many different fields, including chemistry, biology, materials science, and optics [4–6]. In optics, chirality of the matter can be described by the difference in the spectrum when the right-handed circularly polarized (RCP) and left-handed circularly polarized (LCP) waves are incident, i.e., circular dichroism (CD) spectrum [7–9]. Natural chiral materials have limited development in optical device applications due to their weak optical chirality and bulk volume [10]. Artificial nanostructures represented by two-dimensional metasurfaces can control the light field at the subwavelength scale and greatly enhance the interaction between light and matter [11–13]. The interaction with light is governed by the geometry of the meta atoms. With the proposal and development of optical metasurfaces, chiral metasurfaces have attracted widespread attention, and various chiral metasurfaces have been designed and fabricated. Initially, chiral metasurfaces consist of sub-wavelength stacked or helical structures. The chiroptical response of these kinds of metasurfaces is usually caused by the lack of mirror symmetry in the entire three-dimensional (3D) space [14–16]. Although 3D nanostructures offer more freedom in the design process, their fabrication is so difficult that they are limited in large-scale preparation and application. Planar chiral metasurfaces, which only require breaking the in-plane structural symmetry, are the current research focus in this field because of their lower preparation difficulty [17,18]. Traditional planar chiral metasurfaces depend on the local response of the unit cells, leading to low Q factors. Planar chiral metasurfaces with high Q factors can achieve significant chiral field enhancement, which holds great promise for a wide range of applications, such as chiral sensing with high spectral resolution, circularly polarized nonlinear generations, and chiral lasing with low threshold [19]. Due to these advantages, the development of high Q factor chiral metasurfaces is of critical importance.

To overcome this limitation and improve the Q factor while maintaining chiral properties, the introduction of bound states in the continuum (BICs) provides an effective solution [20]. BIC is a localized state spatially and spectrally coexisting with the extended states within the light cone. BICs arise from symmetry mismatch between the BIC and radiated waves, or destructive interference between different radiation channels. They can confine light within the structure for an infinitely long duration, possessing an infinitely large Q factor. BIC can be regarded as a vortex polarization singularity (V point) in momentum space [21–24]. When perturbations are introduced into a BIC system, such as breaking the structural symmetry or altering the incidence angle, the BIC transforms into a quasi-BIC whose Q factor has an inverse square dependence on asymmetry parameters. With special perturbations on BIC metasurfaces, the quasi-BIC at the Γ point in momentum space can be transformed into the C point, which is also called chiral BIC [25]. Chiral BIC metasurfaces can effectively address the low Q factor defects in planar chiral metasurfaces [26–29]. However, designing chiral BIC metasurfaces typically requires extensive parameter scanning, which is not only time-consuming but also computationally expensive. Moreover, since traditional numerical simulation methods are usually optimized for regular structures, this limits the design flexibility and performance of metasurface structures. Recently, the inverse design method has been widely used in the design process of metasurfaces because this new design scheme is more efficient and can break through the limitations of metasurface design and improve the performance of metasurfaces. The inverse design method automatically completes the design of metasurface structures by combining mathematical optimization algorithms such as genetic algorithms, particle swarm algorithms, and neural networks, with electromagnetic simulation [30,31]. The inverse design method can design freeform metasurfaces by specifying the target response rather than the structure of the metasurface itself, which has great design freedom [31]. To date, the application of inverse design methods has been extended to the optimization of various functional metasurfaces. However, the optimization of chiral BIC metasurfaces using inverse design methods has not been studied. Compared with previous chiral metasurface works [32–38], we not only achieved high Q factor resonance under a clean spectral background, but also proposed an automatic inverse design method to replace manual asymmetric parameter scanning. Our method overcomes the limitations of traditional design methods with a global design with a very large parameter space, so that nonlocal effects are considered in the optimizing process, leading to a high Q factor.

In this work, we proposed an inverse design method to optimize chiral BIC metasurfaces. Our inverse design scheme applies an adjoint-based topology optimization method, which is an iterative algorithm that modifies refractive index distribution of a metasurface to maximize its chiroptical response. Through our inverse design approach, a chiral BIC metasurface with 3D intrinsic chirality at the target wavelength is designed and fabricated. The experimental results show that the designed chiral metasurface not only has a strong chiroptical response, but also maintains a high Q factor. Applying our inverse design method on quasi-BIC metasurfaces with different geometric parameters, we can obtain chiral metasurfaces with different response frequencies. Moreover, different fabrication errors are considered in the designing process to improve the robustness of the designed chiral BIC metasurface. We have further extended our method to optimize metasurfaces resonating with elliptic polarizations, which can reflect a specific elliptic polarization state located at the Poincaré sphere while completely transmitting its orthogonal polarization state. The inverse design method we proposed provides an effective solution for the efficient design of chiral BIC metasurfaces.

The flowchart of our design is shown in Fig. 1(a). In our inverse design approach, the quasi-BIC metasurface is taken as the initial structure. This structure has linear eigenpolarizations that break the symmetry of the BIC structure. Then, the quasi-BIC metasurface is optimized using an adjoint-based topological optimization. We define the optimization figure of merit (FoM) as a function of the elements within the metasurface’s Jones matrix. During the iterative evolution, the unit structure of the metasurface continuously updates along the gradient ascent direction of the target FoM, thereby achieving the optimized design of the chiral BIC metasurface. In each iteration, the material gradients at all locations within the meta-atoms are determined by conducting two forward and two adjoint simulations. The detailed algorithm flow is shown in Fig. 1(b) (the detailed topology optimization process is in Appendix A). Once the topological optimization is finished, we obtain a chiral BIC metasurface with specific chiral response under normal incidence. The designed chiral BIC metasurface exhibits total reflection for one circularly polarized light and total transmission for the other, without polarization conversion, as shown in Fig. 1(c).

In the topology optimization inverse design algorithm, a pair of input orthogonal polarization bases is defined, where the incident polarization state and its orthogonal counterpart are defined as |α(ψ,χ)⟩, |β(ψ−90°,−χ)⟩, respectively. The two orthogonal polarization bases can be represented as follows: R(ψ)=[cosψ−sinψsinψcosψ],(1)α=R(ψ)[cosχ−isinχ],(2)β=R(ψ−90°)[cosχisinχ],(3)where ψ and χ represent the orientation angle and ellipticity of the arbitrary polarization states on the Poincaré sphere, respectively; R(ψ) represents the rotation matrix. The output orthogonal polarization bases are also defined as α and β. Thus the Jones matrix J can be obtained by a set of orthogonal polarization bases and their relationships, as follows: J=[rααrαβrβαrββ],(4)where the matrix elements rji (i,j=α,β) represent the conversion coefficients from polarization state i to state j. In order to design a chiral BIC metasurface, the designed metasurface is set to reflect only the LCP light within the target wavelength band, while allowing the RCP light to transmit. Therefore, we use the reflection matrix to characterize the performance of the metasurface, where each element is denoted as Rji=|rji|2. Considering the circular preserving properties of the designed chiral BIC metasurface, the reflection matrix can be written as R=J⊙J*=[RααRαβRβαRββ]=[1000],(5)where ⊙ denotes the Hadamard product; according to reciprocity, the reflectivity Rαβ=Rβα.

Based on the reflection matrix, we can define the FoM in an adjoint-based topological optimization procedure. Consequently, we set the difference between element Rαα and other elements in the reflection matrix as the FoM, which can be written as the following equation: FoM=Rαα−2Rαβ−Rββ.(6)

The variable Rji in the FoM has the gradient of the dielectric constant at each voxel proportional to Re{Eifwd(r) ⋅ Ejadj(r)}, where Eifwd(r) and Ejadj(r) represent the electric field distribution at position r within one unit-cell of the metasurface obtained by forward simulation and adjoint simulation, respectively [9,39]. In the adjoint simulation, the incident polarization states are set as conjugates of |α*⟩ and |β*⟩.

To realize the planar chiral metasurface device with a high Q factor, we use a quasi-BIC metasurface that breaks the in-plane C2 symmetry as the initial pattern for topology optimization. The initial design constructs asymmetric amorphous silicon nanopillars with a height of 320 nm on a quartz substrate. The simulated reflection spectrum is shown in Fig. 2(a). We can find that this structure exhibits a polarization dependent optical response. For x polarized light, the quasi-BIC can be excited at 1307 nm and the full width at half maximum (FWHM) of the quasi-BIC resonance peak is 1.1 nm. Based on this initial design, we further topologically optimized this structure to achieve a specific response to left-circularly polarized light at 1307 nm. During the topology optimization process, it is not necessary to calculate the complete spectrum across the entire target wavelength range in each iteration. Instead, we set the resonant wavelength of the initial quais-BIC structure as the optimization target wavelength. Subsequently, throughout the optimization process, calculations are performed only at this specific wavelength. In Fig. 2(b), we show the changing process of reflection matrix R in the circular basis during the topological optimization. Over the iterations, each optimization parameter gradually approached its target parameter, Specifically, Rll gradually increased and stabilized around 0.78 after 370 iterations, while the Rrr and Rlr(Rrl) continuously declined and finally converged to a minimum. This indicates that the designed metasurface only reflects LCP light at 1307 nm while transmitting RCP light, which aligns very well with the target Jones matrix of the design. The inset of Fig. 2(b) shows that the meta-atom structure pattern is continuously updated with the iterative process in the gradient ascent direction of the predetermined optimization function. During this process, we applied a binarization filter and a Gaussian blur filter to ensure binarization and control the minimum feature size of the optimized structure, ensuring the manufacturability of the optimized structure. In order to further investigate the impact of unit structure changes on polarization response during topology optimization, we also simulated the reflection matrix R in the circular basis of meta-atom structures after 0, 100, 200, 300, and 400 iterations. As shown in Figs. 2(c)–2(g), the chiral response was very weak and the response wavelength is deviated from the design wavelength due to the applied Gaussian blur filter at the beginning of the optimization. As the number of iterations increases, the chiral response of the structure was rapidly enhanced and the resonant wavelength shifted toward the desired wavelength band. When the number of iterations exceeds 200, the reflectivity for the target polarization almost reaches the largest value. At this stage, the reflectivity of the target polarization increases slowly, the FWHM of the resonance decreases gradually, and the other elements of the reflection matrix gradually approach the preset values. Figure 2(g) shows that the Q factor for the optimal structure is 435.

To elucidate the underlying physical mechanisms supported by the designed BIC metasurfaces and their evolution during topological optimization, we calculated the field distribution and multipole decomposition spectrum of the metasurface of the first 100 iterations. By calculating the contributions of various multi-level components to far-field radiation based on induced currents, electromagnetic multipole decomposition in Cartesian coordinates is used. We considered only the electric dipole (ED), magnetic dipole (MD), toroidal dipole (TD), electric quadrupole (EQ), and magnetic quadrupole (MQ), while higher-order terms were neglected due to their minor contributions. These multipole moments can be expressed as follows [40–42]:

electric dipole moment P→=1iω∫j→d3r;(7)electric toroidal dipole moment T→=110c∫[(r→·j→)r→−2r2j→]d3r;(8)magnetic dipole moment M→=12c∫(r→×j→)d3r;(9)magnetic toroidal dipole moment G→=k220∫r→2(r→×j→)d3r;(10)electric quadrupole moment Qab(e)=12iω∫[rajb+rbja−23(r→·j→)δab]d3r;(11)magnetic quadrupole moment Qab(m)=13c∫{(r→×j→)arb+[(r→×j→)bra]}d3r,(12)where ω denotes the angular frequency, r→ is the position vector, k is the wavevector, c is the speed of light, and a, b=x, y, z. Here, we ignore the small contribution of higher-order multipoles. The scattered powers of these multipole moments are calculated from Ip=2ω43c3|P→|2,IM=2ω43c3|M→|2,IT=2ω63c5|T→|2,IG=2ω63c5|G→|2,IQ(e)=ω65c5|Qab(e)|2,IQ(m)=ω640c5|Qab(m)|2.(13)

The total scattered power of the multipole moments can be given as Itotal=IP+IM+4ω53c4(P→·T→)+IT+4ω53c4(M→·G→)+IG+IQ(e)+IQ(m),(14)where 4ω53c4(P→·T→) and 4ω53c4(M→·G→) are the interference terms of the ED and ETD, MD and MTD. The multipole decomposition spectrum of the optimized structure during the optimization process is shown in Figs. 3(a)–3(e). The corresponding electric field and magnetic field distributions are shown in Figs. 3(f)–3(j) and Figs. 3(k)–3(o), respectively, where white arrows represent the magnetic vector and black arrows represent the induced current. The initial structure supports a quasi-bound state in the continuum (quasi-BIC), with the TD dominating at the resonant wavelength. As shown in Figs. 3(f) and 3(k), the current vector distribution in the xz plane forms a toroidal pattern, generating an MD oriented along the y direction. Two oppositely directed MDs connected end-to-end in the xy plane further produce a TD along the z direction. To ensure smooth structural variations, we incorporated Gaussian blurring into the topological optimization process. As iterations progressed, the geometry of the meta-surface unit cells evolved, leading to changes in the electromagnetic field distributions. As illustrated in Figs. 3(g) and 3(i), magnetic field distribution along x direction also occurs in addition to the induced current loop in the xz plane, which leads to the enhancement of MD in Fig. 3(a). As shown in Figs. 3(m)–3(o), the magnetic field distribution forms a loop in the xy plane within 20–100 iterations. At the same time, the induced current changes, as shown in Figs. 3(h)–3(j). Throughout the optimization process, the TD remained dominant at the resonant wavelength, shown in Figs. 3(a)–3(e). As the number of iterations increased, the contribution of the TD gradually intensified, accompanied by a blue shift of the resonance peak toward the target wavelength. The optimized structure exhibits significantly enhanced TD contributions at the resonant wavelength, along with a blue shift of the resonance peak (the results of 0–400 iterations are in Appendix B).

Finally, we conducted experiments to verify the feasibility of the chiral BIC metasurface design. We fabricated three metasurfaces, each with a size of 200μm×200μm, realizing the functions of BIC, linearly polarized chiral BIC, and circularly polarized chiral BIC at the target band. The three metasurface samples were fabricated on a 320 nm thick amorphous silicon layer transferred onto a 1 mm thick glass substrate. The patterns were defined using electron beam lithography (EBL) and reactive ion etching (RIE) processes. In the near-infrared wavelength range, the imaginary parts of the refractive indices of silicon and glass are very small, and the absorption losses of the materials can be neglected (see Appendix C for details). Therefore, we conducted verification experiments using a transmission optical path to minimize experimental errors, as shown in Fig. 4(a). The light from the broadband laser source (NKT Photonics SuperK COMPACT) was focused onto the fabricated metasurface using a lens and then collected using an optical spectrum analyzer (Yokogawa AQ6370C) through a fiber coupler. The measured transmission spectra of the fabricated three BIC metasurfaces are shown in Figs. 4(b) and 4(c), and the corresponding scanning electron microscope (SEM) images are shown in the insets.

The transmission spectra obtained from the manufactured metasurfaces are essentially consistent with the calculated results presented in Figs. 2(a) and 2(e), exhibiting the polarization characteristics of BIC metasurfaces. Figure 4(b) shows that the linearly polarized quasi-BIC is characterized by a resonant wavelength of 1300 nm, FWHM = 9.7 nm, and Tx=0.23. Figure 4(c) shows that the circularly polarized chiral BIC is characterized by a resonant wavelength of 1298 nm, FWHM = 13 nm, and Tll=0.39. The resonant peak of the metasurface shifted from 1307 nm to 1298 nm, and there is a deviation in the overall transmittance and Q factor compared to the simulated design. To investigate the optical characteristics of the topologically optimized metasurface, we analyzed the CD under normal incidence of RCP/LCP waves. We define a CD as the transmittance difference under LCP and RCP incidence: CD=(Tll+Trl)−(Trr+Tlr)(Tll+Trl)+(Trr+Tlr).(15)

Figure 4(d) shows the simulated and measured CD spectra. The minimum value in the CD simulation is −0.76, while the minimum value in the experiment is −0.32. The reduction in spectral intensity and the broadening of the FWHM between the simulation and experimental results are primarily attributed to fabrication imperfections.

To demonstrate the versatility of our proposed method, we also utilized the same process to fabricate metasurfaces with resonance peaks at other wavelengths. We optimized the chiral BIC metasurface around 1225 nm. Figure 5(a) represents the simulated spectral curves of the optimized metasurface. The corresponding experimental spectra are shown in Fig. 5(b), which show the chiral responses at 1245 nm. As shown in Fig. 5(c), the experimentally measured CD values exhibit a broader bandwidth at the responsive wavelength compared to the simulated CD values. We believe that the discrepancies between the simulation and experimental results are primarily attributed to manufacturing defects and errors (see Appendix D for details). These results demonstrate the effectiveness of our design methodology and confirm the manufacturability of the metasurface devices designed using this approach.

In addition to designing chiral metasurfaces, we have further extended our method to optimize metasurfaces resonating with specific elliptic polarizations. This metasurface allows preferential reflection of the polarized state α(ψ,χ) and full transmission of its orthogonal state β(ψ−90°,−χ). First, we design metasurfaces with polarization along the latitude lines of the Poincaré sphere through topological optimization. Then, by physically rotating the metasurface, we can achieve coverage of longitudes on the Poincaré sphere, which corresponds to polarization states with the same ellipticity χ but different orientation angles ψ. By combining these two steps, the polarization of the optimized metasurface can be almost anywhere on the Poincaré sphere. We used the aforementioned meta-atom structure exhibiting a quasi-BIC response at 1307 nm as the initial design for topology optimization and designed four sets of quasi-BIC metasurfaces with ψ=0° and 2χ=±π/3, ±π/6. The positions of these polarization states on the Poincaré sphere are shown in Fig. 6(a). In the optimization for elliptic polarization, we set the output orthogonal polarization bases as α and β, the same as input. Figures 6(b)–6(e) display the meta-atom structures and their performance, which have a reflection rate of over 0.8 for the target polarized light and a response below 0.1 for the orthogonal polarization.

In summary, we have proposed an inverse design approach based on topological optimization to enhance chiral BIC metasurfaces. Using the proposed inverse design method, we have theoretically and experimentally verified circularly polarized chiral BIC metasurfaces with circular-polarization-preserving components under normal incidence. The designed chiral BIC metasurface can reflect the left-handed circularly polarized component while transmitting its orthogonal component without circular polarization conversion. The optical response of the fabricated chiral BIC metasurface is consistent with the design. The fabricated chiral BIC metasurface shows a resonance of Q∼100 and a minimum CD of −0.32. Furthermore, to enhance the versatility of the algorithm, we have simulated and verified its applicability to the design of quasi-BIC metasurfaces for specific elliptic polarization. Additionally, the algorithm can be further extended to design multi-wavelength chiral BIC devices. This method is not only theoretically innovative but also shows broad potential in many applications, especially in the fields of chiral imaging, chiral sensing, high-selectivity optical filters, and chiral metasurface holographic displays.

We employ inverse design topology optimization to engineer metasurfaces with specific polarization behavior, enabling the generation of freeform structure. The numerical simulations are performed using the RETICOLO software package based on the rigorous coupled-wave analysis (RCWA) algorithm [43]. Our inverse design approach employs gradient ascent to maximize the figure of merit (FoM). The reflection matrix for specific elliptic polarizations is obtained by {α,β}−1Jlin{α,β}.

For this fabrication process, we deposited a 320 nm thick amorphous silicon film on a 1 mm double-polished glass substrate using an ICP-PECVD system (SENTECH SI500D). The patterning was carried out on a positive photoresist (PMMA 495 A4) using an electron beam lithography (EBL) system (Raith e-lineplus). After development and fixing, the pattern was transferred onto a 30 nm thick chromium layer via electron beam evaporation. Chromium was then used as a hard mask for etching the amorphous silicon layer through a reactive ion etching system (SENTECH SI591). Lastly, residual chromium was removed using an ammonium cerium nitrate solution.