Metasurfaces, composed of localized sub-wavelength structures, enable arbitrary control of electromagnetic wavefronts. Initially, discrete structures receive great attention for realizing phase¹, amplitude^2,3, and polarization manipulation^4-6, and through the precise control of electromagnetic waves, leading to many useful applications, such as focusing^7-9, polarization conversion^10-12, filtering¹³, holographic imaging¹⁴, encryption^15-17, and information manipulation^{18, 19}. Based on the diffraction theory, the diffraction efficiency limit is determined by the discrete level N of quantized phase distribution. For N values 2, 4, 8, and 16, the corresponding diffraction efficiency limits are 40.5%, 81.1%, 95.0%, and 98.7%²⁰, respectively. To enhance the diffraction efficiency, it is possible to increase the value of N. However, it will lead to extra electromagnetic coupling between adjacent meta-atoms and manufacturing difficulties.

In 2015, a quasi-continuous structure, namely a catenary structure, was proposed to overcome the challenges associated with discrete geometric phases²¹. The optical devices that impart a continuous phase to the incident light are known as “catenary diffractive waveplates”²². Since then, catenary optics has emerged as a new research direction^23-31. Due to the continuous phase control ability of the catenary-like structure, it has an advantage of higher efficiency for broadband modulation and is widely applied in high-performance devices, such as efficient beam acceleration generation²⁹, tunable circularly polarized beam splitting³², and polarization-controlled unidirectional excitation of surface plasmon polaritons³⁰. Compared with metallic catenary structures, dielectric catenary structures exhibit lower material losses²⁵, resulting in a higher diffraction efficiency. The diffraction efficiency limit of metallic catenary structures is only 25%, while dielectric catenary structures surpass this limit^{23, 27}. However, the performance of nanostructures in dielectric materials is sensitive to their dimensions²⁵. Especially for equal-width catenary structures, the inconsistent filling ratio of local structures may lead to an undesired parasitic propagation phase, leading to a reduction in diffraction efficiency²⁷. In previous research, a topology optimization method was proposed to further optimize geometrical structures into quasi-continuous structures³³. The results demonstrate that enhancing the degree of freedom in structural design can significantly improve the devices’ performance. Nevertheless, high-degree-of-freedom structures require a large amount of data for characterization, leading to challenges in the design of large-scale devices. Unlike topology optimization, adjoint-based shape optimization takes the structures’ boundaries as optimization parameters, enabling the rapid design of large-scale devices due to its straightforward parameterization characteristics. However, adjoint-based shape optimization is currently only used to optimize the propagation phase devices^34-36.

In this work, we proposed a novel approach named field-driven optimization (FDO) method to improve the performance of catenary optical devices. Such a method utilizes real boundary conditions to accurately calculate the electromagnetic field in each iteration, facilitating precise control of both geometric and propagation phases of the structures. By iterating the widths, the local optimal structure can be achieved with only 10 iterations of optimization. More importantly, the FDO method has the advantage of low parameterization cost, which enables a small number of parameters to characterize the structure. Above all, the establishment of a structure library containing catenary structures with varying periods enables the rapid design of centimeter-scale devices using a simple boundary interpolation method. This approach effectively overcomes the limitation of computer memory in simulation calculations. In addition, we have also achieved efficient phase control for incident light with arbitrary polarization states using the proposed FDO method. Our work proposes the FDO method to optimize the catenary optical devices for the first time, which may further promote the development of catenary optics.

To achieve the performance limit of catenary optical devices, it is imperative to precisely modulate both the geometric and propagation phases. The geometric phases originate from the photonic spin-orbit interactions (SOI), which are dependent on the orientation angle θ^37-39, while the propagation phases are size-dependent⁴⁰. Here, the fabrication of catenary structures is achieved by the integration of forward design and adjoint optimizations. As illustrated in Fig. 1, the trajectory (solid yellow line) and the width parameter w of these structures completely define the design. By leveraging the Jones matrix theory, the trajectory of the structures is obtained through forward design. Then, the trajectory is used in conjunction with initial widths, serving as the initial structures for subsequent optimization using the FDO method.

Firstly, the trajectory of the structures can be integrated through the computation of the orientation angle θ of each point along the trajectory. Ideally, the orientation angle of each point can be calculated using the ideal Jones matrix J_i, which can be written as:

Ji=eiφ1q1*q1†+eiφ2q2*q2†,（1）

where q₁=[cosχ sinχe^iδ]^T and q₂=[sinχ −cosχe^iδ]^T represent two arbitrary orthogonal states of polarization in the linear polarization basis, and χ and δ determine the respective polarization states. φ₁ and φ₂ are ideal phases for q₁ and q₂, respectively. The superscript * is complex conjugation, T denotes transpose operation, and † indicates the Hermitian conjugate. The eigenvectors of J_i determine the required orientation angle θ of each point. The trajectory of the actual structure is integrated by y=∫tanθdx. Then, as illustrated in Fig. 1, the width w is defined as the distance along the trajectory’s normal direction. Thus, the coordinates of the boundaries are defined as:

xu,l=x∓wu,lsinθyu,l=y±wu,lcosθ,（2）

where (x_u,l, y_u,l) indicates the coordinate of the upper and lower boundary, (x, y) is the corresponding coordinate of the trajectory, and w_u,l is the width between the upper, lower boundary and the trajectory along the normal direction. In previous research, the width of a single catenary structure is often a constant, which will introduce an additional propagation phase gradient, limiting the efficiencies of the devices. To achieve near-perfect modulation of both geometric and propagation phases of optical devices, the FDO method is proposed to correct the propagation phases.

Based on the adjoint theory, the FDO method is an iterative algorithm that modifies the widths of the continuous structures to maximize the figure of merit (FoM)³⁴. Compare to traditional optimization methods such as genetic algorithms^{41, 42}, evolutionary algorithms^{43, 44}, and particle swarm optimization algorithms^{45, 46}, which rely on random exploration, the proposed approach can obtain gradients of all variables relative to the FoM using only two simulations, namely the forward and adjoint simulations^{33, 35, 36, 47}. Figure 1 shows an FDO optimization example of a one-dimensional (1D) metalens composed of catenary-like structures. The optimization process starts by defining the FoM. The goal of the designed device is to transform the normal incident light field into the ideal field distributionE_ideal. Thus, the FoM could be defined as the fraction of the output field in the ideal field distribution³³:

FoM=|〈Eideal|Ef〉|2,（3）

where E^f is the field realized by the current design in the forward simulation. The subsequent imperative step requires conducting another simulation, namely adjoint simulation, to get the adjoint field E^A. Different from adjoint-based shape optimization, which optimizes the length and width of the whole structure, the proposed FDO method optimizes the width of each point on the trajectory of the structure iteratively. The gradient of the widths of each point on the trajectory could be set as^{36, 48}: (The details are shown in Fig. S1, Supplementary information.)

δwu,l=Re{(ϵ2−ϵ1)[Ex′A(xu,l,yu,l)Ex′f(xu,l,yu,l)+Dy′A(xu,l,yu,l)Dy′f(xu,l,yu,l)ϵ1ϵ2]}[Ex′Ey′]T=R(θ)[ExEy]T,（4）

where ε₁ and ε₂ denote the dielectric constants of the environment and the structures, respectively. E_x′,y′ and D_x′,y′ represent the components of the electric field and electric displacement vector along the x′, y′ direction, which is tangent and normal to the trajectory of the structure as shown in Fig. 1. E_x and E_y are the x, y component of the electric field obtained by the two simulations, and R(θ) denotes the coordinate transformation matrix.

It’s worth noting that our method is applicable to incidence with arbitrary polarizations. The FDO theory is polarization-independent. During the process of initial structures design, the incident polarization state can be altered by changing the q₁ and q₂.

To demonstrate the effectiveness of the proposed method, the performance of the catenary array for circularly polarized (CP, χ=45°, δ=90°) incidence is optimized. The trajectory of the catenary remains unchanged to maintain the geometric phase gradient of the structures, while the widths of the catenary are optimized iteratively to correct the propagation phase gradient. In addition, to extend catenary optical devices for incidence with arbitrary polarization, the periodic catenary-like arrays are also designed and optimized for elliptically polarized (EP, χ=50°, δ=60°) and linearly polarized (LP, χ=0°, δ=0°) incidences, respectively. Ideally, these three different optical devices are designed to deflect the incident CP, EP, and LP lights to the +1^st order at the central wavelength of 10.6 μm, using pure geometric phases, co-regulation of geometric and propagation phases, and pure propagation phases, respectively. Assuming the phase gradients attached to a pair of orthogonally polarized light by the device are opposite, the ideal phases are φ1,2=±2πsinαλx, where α is the diffraction angle, and λ is the working wavelength. For simplicity, the widths associated with the points along the trajectories of the initial structures are either constant (for CP and EP incidences) or linearly proportional to the x-coordinates of the trajectory (for LP incidence). According to Eqs. (1) and (2), the initial structures of the optical devices can be designed for arbitrary polarized light illumination. The adjoint source is set as an ideally polarized light with a phase distribution of −2πsinαλx+φ0, and then conjugated, where φ₀ is a phase constant that is updated based on the phase distribution of the forward simulation field in each iteration. The electromagnetic simulations are performed using the finite difference time domain (FDTD) method.

Figure 2(a) depicts the boundary variations of the three devices with a diffraction angle of α=45° for CP, EP, and LP incidences, respectively. The dimensions of the initial and shape-optimized structures are shown in Fig. S2, Supplementary information. The materials of the substrate and the structures are set as silicon (Si) with an index of n=3.37. To avoid the intersection between adjacent structures, the catenary-like structures for CP and EP lights are truncated on both sides. The optimization process is demonstrated by the evolution of the diffraction efficiencies of the optical devices, as illustrated in Fig. 2(b). The diffraction efficiency is defined as the power ratio between the transmitted light at the specific diffraction order and the total transmitted light. As shown in Fig. 2(b), the efficiency of the catenary optical device with α=45° increases from 96.4% to 99.2% due to the modulation of the propagation phases, achieving near-perfect beam wavefront modulation. In fact, with 10 iterations, the efficiency of the catenary optical device is already very close to the optimal value, reaching 99.5% of the optimum. For the catenary-like optical devices designed for EP and LP incidences, the diffraction efficiencies eventually improve from 87.9% and 92.0% to 94.2% and 98.6%, respectively. The bottom one of Fig. 2(b) shows the evolution of the diffraction efficiencies of the devices with α=75°. After 30 iterations, the diffraction efficiencies of the three devices all exceed 90%. Among them, the efficiency of the catenary device reaches up to 94.5%. Although the devices are primarily optimized at the wavelength of 10.6 μm, a significant enhancement in performance in observed across the broadband range of 8–13 μm. The details are shown in Fig. S3, Supplementary information.

Unlike periodic optical devices, the coupling effect between adjacent structures also affects the performances of devices with non-periodic distribution. Here, two non-periodic 1D metalenses are designed for CP and LP incidences using the FDO method. The ideal phase distributions of the metalenses are defined as: −2πλ(x2+f2−f), where f=850 μm denotes the focal length of the metalenses. The diameters of the metalenses are D_lens=1000 μm, and the working wavelength is 10.6 μm. For CP incidence, the geometric phases are used to modulate the wavefront of the incident light. To achieve the ideal geometric phases, the equal-width catenary-like structures with width w=1.1 μm are used as the initial structures. Each catenary-like structure (including the rectangle regions) is utilized to realize the geometric phases coverage from 0–2π, and their trajectories are integrated in accordance with the ideal geometric phases. Different from the devices that use geometric phase modulation with CP light incidence, optical devices employ propagation phases to modulate LP light. As a result, the initial structures of the device feature a trapezoidal distribution.

The efficiencies of the metalenses for CP and LP incidences are depicted in Fig. 3(a) and 3(d), respectively. After 10 iterations, the efficiency trends of the two metalenses have substantially converged. Finally, at the 30th iteration, the focusing efficiencies of the two metalenses improve from 77.39% and 76.95% to 91.71% and 91.28%, while the absolute efficiencies enhance from 62.18% and 58.82% to 78.55% and 72.66%, respectively. Here, the focusing efficiency is defined as the ratio of the power within the range of three times the full-width half-maximum (FWHM) and the energy of the transmitted light, divided by the ideal one. The ideal light intensity distribution can be obtained using the vectorial angular spectrum theory of a polarized wave carrying the ideal phases⁴⁹. The FWHM of the ideal intensity distribution on the focal plane is 9.35 μm for CP incidence, while it is 9.77 μm for LP incidence. Absolute efficiency defines as the product of the focusing efficiency and the transmittance. Figure 3(b) and 3(e) show the normalized light intensity distributions of the initial and shape-optimized metalenses on the focal plane. The insets in Fig. 3(b) and 3(e) demonstrate that the cross-sections of the light intensity profiles of the shape-optimized metalenses exhibit less stray light, providing that the shape-optimized metalenses can effectively converge the incident light onto the focal spot. Compared to the initial metalenses, the light intensity distributions of the shape-optimized ones are closer to those of the ideal Airy spots. The maximum light intensities at the focal spot are enhanced by 23% and 45% for CP and LP incidences, respectively. Assuming that the metalenses are placed on the XoY plane and the incident light is propagating along the +z direction, Fig. 3(c) and 3(f) show the distributions of light intensity on the XoZ plane, simulated by FDTD solutions.

We proposed a novel boundary interpolation method to overcome the limitation of computer memory on the diameter of the optical devices that need to be optimized. This method enables the fast design of large-scale complex devices with high performance. The ideal phases of the optical device with a particular polarized incidence are determined by the propagation phases and geometric phases. Among them, the geometric phases can be strictly obtained by Eq. (1), while the propagation phases correspond to the geometric phases. Due to the parametrizable of the shape-optimized structures, the relationship between the geometric phases and propagation phases (i.e., the widths of the structures) can be easily established. To design large-scale optical devices, similar to the design approach mentioned above in Section Theory, the geometric phases are first obtained according to the ideal phase distribution. However, the widths of each structure are no longer obtained using the FDO method but are predicted according to the relationship between the geometric phases and widths of the structures in the library. The widths of the large-scale devices can be obtained without requiring the simulation of the electromagnetic field, allowing the fast inverse design.

To enable the rapid design of optical devices illuminated by CP light with a center wavelength of 10.6 μm, a catenary library is established via the FDO method. Figure 4(a) shows the shape-optimized catenary structures in the library. The local width of the predicted structure is obtained by interpolation based on relevant parameters, including the diffraction angles, orientation angles, and their corresponding local widths, of the structures in the library. Therefore, the more the optimized structures in the catenary library, the closer the local width distribution of predicted one will be to those obtained using FDO method, resulting in a better performance of the predicted structures. To substantiate the method, our library comprises 9 shape-optimized catenary structures with different diffraction angles (α=5°, and 10° to 80°, with an interval of 10°). These 9 structures could cover diffraction angles ranging from −90°~90° (the negative diffraction angles can be realized by inverting the corresponding catenary structures), facilitating a rapid design of devices with arbitrary diffraction angles. The number of shape-optimized structures in the library can be further increased to achieve devices with enhanced performance. The initial structures are set as equal-width catenary structures with an initial width of 1.1 μm. The complex conjugate of the adjoint source is set as an ideally polarized plane wave with phase distribution −2πsinαλx+φ0. Unlike the iteration process described in Section Periodic catenary-like devices for arbitrary polarization of light, here φ₀ is a fixed constant in all iterations to ensure the propagation phases attached by each catenary structure are close. To facilitate the subsequent design of high-performance optical devices using the interpolation method, we first established the relationship between the arbitrary orientation angle and the width of the local structures through interpolation. As shown in Fig. 4(a), the continuous part of the catenary structures is indicated by the inner dashed box, while the rectangles part of the catenary structures is represented by the outer box. These discrete portions of the structures may result in abrupt changes in the widths of a single catenary structure.

Theoretically, the library can be used to design arbitrary geometrical phase optical devices that operate at a wavelength of 10.6 μm. Based on the library, firstly, the geometric phases are calculated by Eq. (1). Then, the trajectories of the structures are designed based on the geometric phases. The widths of the structures are obtained using a two-dimensional (2D) interpolation method based on the horizontal lengths and geometric phase distributions of the desired structures and those of the structures available in the library. Accordingly, catenary optical devices with high performances are designed. Based on the library, 7 catenary arrays with different diffraction angles (from 15° to 75°, with an interval of 10°) are predicted. As shown in Fig. 4(b) and 4(c), the efficiencies of the predicted catenary devices are close to those of the devices designed using the FDO method, exhibiting an average difference of less than 2% in both diffraction and absolute efficiencies. Besides, compared to the initial devices with equal-width catenary structures, the efficiencies of the predicted ones are significantly improved, with all diffraction efficiencies surpassing 94%, and absolute efficiencies increasing by ~8%. The absolute efficiency of the deflector is defined as the product of the diffraction efficiency and transmittance. As shown in Fig. 4(a), there is a sharp change in the local width of the catenary structure near the orientation angle of 0°, which corresponds to a diffraction angle of 30°. Such a sharp change may be a result of the FDO process getting trapped in a local optimum, and lead to inaccuracies in predicting catenary structures with a diffraction angle of 35°. Increasing the number of structures in the catenary library might potentially address this issue. In addition to the devices with catenary structures, the catenary library shows significant advantages in designing optical devices with catenary-like structures as well. Figure 4(d) illustrates the comparison between the focusing efficiencies of the 1D metalenses (D_lens=1000 μm) with equal-width catenary-like structures and those of the metalenses designed using the library with numerical apertures (NA) ranging from 0.1 to 0.7. The results demonstrate that the performance of the predicted metalenses improved by ~15% compared to that of the equal-width ones. In addition to the fast design, our approach splits the optimization of large-scale devices into the optimization of several periodic arrays, which renders the design problem computationally tractable.

For further verification, we designed the 1D catenary-like metalenses with diameters of 6000 μm using the library and measured their performances. The experimental setup is shown in Fig. S4, Supplementary information. Figure 5(a) shows the local scanning electron microscope images of the metalens with NA=0.3. The experimental focusing efficiencies and zero-order efficiencies are shown in Fig. 5(b) and 5(c). Here, the experimental focusing efficiency (zero-order efficiency) is defined as the energy ratio between the focused light (zero-order light) and the total transmitted light. The “zero-order light” refers to the diffracted light that propagates along the optical axis of the incident light. To measure the experimental focusing efficiency and zero-order efficiency, the power detector is placed at the focal plane of the metalens at distances of d ≈ 0.5 cm and 25 cm, respectively. The details are shown in Fig. S4, Supplementary information. The focusing efficiencies of the catenary-like metalenses designed based on the library have shown significant improvement of ~15% compared to the equal-width ones, whereas the zero-order efficiencies have decreased by ~50%. To further visibly compare the performances of the metalenses, we experimentally evaluate the performance of the catenary-like metalenses with NA=0.3, 0.4, and 0.5. The measurement setup is shown in Fig. S4, Supplementary information. A charge coupled device (CCD, pixel size: 17 μm × 17 μm) camera is used to detect the light intensity distribution of the focal plane. Based on the vectorial angular spectrum theory, the theoretical FWHM of the focal spot for the metalenses with NA=0.3, 0.4, and 0.5 is approximately 15.70 μm, 11.81 μm, and 9.48 μm, respectively. The experimental light intensity distribution at the focal spot is under sampled. Thus, we conducted a comparison of the background noise at the focal plane of the metalenses overexposure. Figure 5(d–f) show the normalized intensity distributions of the metalenses with different NA on the focal plane. During the measurement process, the power of the incident light and the parameter settings of the CCD are kept unchanged. Figure 5(d–f) show that the predicted catenary-like metalenses significantly reduce the stray light on the focal plane by ~40% in comparison to metalenses composed of equal-width catenary structures. Noting that it takes only 2.0 s to design such a metalens via the library (processor: Inter(R) 3.40 GHz; memory: 8G). Directly optimizing the full mode of such a device is extremely difficult. Our computer faces the problem of memory shortage when working on such device using FDTD solutions. (CPU: Inter(R) 2.30 GHz; memory: 256 G; mesh: 2 μm). More intuitively, our method can achieve the design of a metalens with a diameter of ~2000λ in 8 s, which reduces the optimization time by several orders of magnitude compared to that of the advanced adjoint-based algorithms available today⁵⁰.

To achieve the performance limit of catenary optical devices, precise phase modulation is employed. Geometric phases are strictly calculated based on the optical Jones matrix, while the propagation phases are corrected using the FDO method. The catenary structures can be optimized to achieve perfect wavefront modulation, with a resulting diffraction efficiency near 100% in 10 iterations. Besides, our work expands the catenary optical devices to incidence with arbitrary polarization. Although our approach currently only demonstrates optimization for 1D metalenses, 2D metalenses can be designed and optimized by modifying the integral path of the catenary-like structures and adjusting the phase distribution of the metalenses, while the entire optimization progress remains the same. To enable the rapid design of large-scale high-performance optical devices, we establish a comprehensive library of parametrizable shape-optimized structures, which greatly facilitates the rapid design of high-performance optical devices over a large scale. Simulation and experimental results show that devices designed by the library of only nine curved structures significantly improve performance by ~15% compared to devices with equal-width structures. The addition of more structures to the library would increase the available options and enable finer adjustments, resulting in more precise interpolations and better-performing devices. In addition, although the library is currently limited to design devices for CP light only, the method has the potential to establish a library for arbitrarily polarized incidence. In the future, it may be possible to achieve rapid design of devices with arbitrary polarized light incidence by establishing a library that contains the information about the incident polarization state. Besides, combining this method with deep learning may further improve the efficiencies of the devices.