{'fraction_non_alphanumeric': 0.027109621451104102, 'fraction_numerical': 0.01015378548895899, 'mean_word_length': 4.76228344220392, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Metasurface Hybrid Lens Design using a Semi-Analytical Model

We propose a new method for integrating metasurfaces in optical design using semi-analytical modelling of dielectric nanostructures. The latter computes the output phase of an electric field incident on the metasurface, allowing their use with ray-tracing software. This tool provides a method to use metasurfaces in optical systems while using built-in optimization processes to avoid time-consuming computation. To demonstrate the applicability and versatility of our method, we present variations of a triplet composed of refractive elements and a metasurface. For each of the systems, similar optical performances are achieved. Our unique and innovative approach to joining metasurfaces and ray tracing has the potential to promote the design of innovative systems by exploiting the richness of metasurfaces and the functionality of conventional lens design software.


Introduction
Research on metasurfaces has been motivated by their broad potential for manipulating light. Exotic applications exploiting their increased control on the wavefront, their polarization sensitivity or their unique chromatic dispersion were proposed in the past years [4,1,3,17,12]. Many of these technologies were demonstrated by wave optics computation, sometimes followed by experimental work. Indeed, high nonlinear responses [16], multifunctional metalens [23,11] or non-Hermitian metasurfaces [9] can be found in the literature, to cite a few.
The utility of metasurfaces has been praised many times; some even say they have the potential to replace traditional optics [20,26,14]. However, their integration in conventional lens design remains a little-explored avenue. Considering metasurfaces in a complete optical system is inconceivable if common time-consuming numerical methods are used. The macroscopic scale of elements and the large number of iterations required to design a suitable system [24] would require too much computing power. Therefore, a method to efficiently represent metasurfaces in a ray-tracing process is needed. Such a method could originate from modelling the optical properties of meta-atoms.
Thus, we propose an approach for using metasurfaces in ray-tracing applications thanks to the efficient characterization of their optical properties provided by a model. Section 2 describes the ray-tracing equations adapted for propagation through metasurfaces. Such equations will later be implemented in ray-tracing software like Zemax ® . Then, in section 3, a semi-analytical model computing the optical properties of a metasurface composed of dielectric nanocylinders is given. This modelling is a continuation of our previous work, where rectangular structures were studied [2,6]. Results obtained by this model are then compared to FDTD simulations in section 4. Finally, we present optical designs using different refractive lenses and metasurfaces targeting the same application to display the method's versatility in section 5. arXiv:2209.06041v1 [physics.optics] 13 Sep 2022 2 Ray-tracing through a metasurface As the ray-tracing method is based upon calculating the angle of refraction at each optical surface, it is natural to start our implementation with the generalized law of refraction for metasurfaces proposed by Yu et al. [28]. Here is the case for the general refraction where the metalens implements a phase ϕ on the wavefront, and n i , θ i and n t , θ t are the refractive index and the incidence angle before and after the metasurface.
In optical design software such as Zemax ® or CodeV ® , rays are represented by a three-dimensional vector R = xx + zẑ + zẑ. The orientation of such a ray is defined by its direction cosines l, m and n. l = cos u x , m = cos u y , n = cos u z , where u x , u y et u z are the angles formed by the vector R and the axesx,ŷ andẑ. Considering metasurfaces are flat components placed perpendicular to the optical axis in most cases, we can establish that θ i = 90 • − u =⇒ sin θ i = cos u. Thus, by rearranging Eq. 1 with n t = n i = 1 using the direction cosines, the following set of equations for the coordinates (l , m , n ) of the refracted beam can be obtained The normalization condition l 2 + m 2 + n 2 = 1 was used for the n -coordinate. If one wishes to perform ray tracing through a metasurface, one must know the phase shift induced for a given incident beam 3 Semi-analytical model As stated previously, obtaining the phase shift implemented by the metasurface via conventional full-wave simulations is inconceivable for lens design purposes considering the standard optimization process. Therefore, a model that can efficiently output the optical properties used for ray-tracing is advisable. Therefore, this section concerns the development of an expression for the phase shift of cylindrical nanostructures, an essential parameter for performing the ray tracing as previously determined.
Related to our earlier work about nanofins [2], the guided propagation through a nanocylindrical waveguide can be described as a phase accumulation of the electrical field given by the optical path of an equivalent wave in a medium of length H as a first-order approximation.
where τ j and β j are respectively the complex amplitude and the propagation constant of the waveguide mode j. For the needs of our application, we consider the fundamental mode, leaving only to obtain the complex amplitude and the propagation constant of the latter.

Propagation constant
The propagation constant β of a single nanocylindrical structure can be obtained by solving the following dispersion relation [27,5]: with p 2 1 = ω 2 µ 0 1 − β 2 , q 2 2 = β 2 − ω 2 µ 0 2 . ω, µ 0 , 1 and 2 are the common physical constants. J n and K n are the Bessel function of the first kind and the modified Bessel function of the second kind order n respectfully with prime as the derivative of those function.
From the propagation constant of an isolated structure, it is necessary to add a correction to consider the coupling effects, knowing that a metasurface is composed of an array of equidistant nanostructures. To do so, Coupled Mode Theory (CMT) is often used because of its simplicity, but multiple approximations should be met for its computation to be accurate [10]. To consider the periodicity, we will apply a correction to the propagation constant following the Minot method [19]. Using the slowly-varying envelope and assuming a coupling C which is much smaller than the propagation constant, a long-established approximation of the effective propagation constant of an array of waveguides is expressed as where S the distance between waveguides, k x is the transverse Floquet-Bloch waves and C the coupling. Here, we will assume that k x = p+1 N +1 π for the p structure within an N -size array [19]. The coupling C can be computed following with 00 and 001 expressed as where M p correspond to the transverse field distribution of the p structure of the array. For our implementation, N = 2 and p = 0 are used, corresponding to the correction of the closest neighbour. For more details on the calculations applied for cylinder structures, see appendix A.

Complex amplitude
Using the Snell-Descartes law, one obtains the transmission angles at the first and second sides of the nanocylinder. Because of the high contrast of the refractive indices often used in dielectric metasurfaces [3,13], multiple scattering events might occur. To take this phenomenon into account, we consider the Fabry-Perot effect of the nanostructure by introducing a simple geometric sum formulated as where t 1,2 , r 1,2 correspond to the transmission and reflection coefficients at the first or second interface of the nanofin, and β is the propagation constant inside of the nanocylinder of height H.
Finally, from Eq. 6, the output phase-shift ϕ produced by a nanorod of radius a considering only the fundamental mode is given by

Meta-atoms properties
The semi-analytical model allows for computing the phase shift introduced by those nanostructures under various illuminations. Fig. 1a presents the phase shift ϕ for 600 nm high TiO 2 nanocylinders on a SiO 2 substrate for wavelengths λ from 568 nm to 628 nm. The structures at a varying radius between 45 and 110 nm are located 350 nm apart. These parameters were used since they correspond to metasurfaces often represented in the literature [18,7,22]. As one can observe, the phase can be accurately modelled for various wavelengths when compared to FDTD simulations using Lumerical software [25]. We can also extend our analysis to illuminations at various incidence angles, as shown in Fig.  1b where the incident beam of light comes at a 20 • field of view. To assess the magnitude of the error in our predictions relative to numerical simulations and the impact this would have in a lens design context, we refer to the work of Patoux et al. [21]. Through their research regarding nanofabrication errors and their effect on metasurfaces performances, they concluded that the fabrication error must be limited to 10 nm for a metalens to maintain acceptable quality. We have identified that the deviations observed in Fig. 1 are equivalent to a variation of 1-4 nm in radius on average which is well below the specified limit of 10 nm. Furthermore, a numerical study on the variation of the Strehl ratio showed that errors similar to those shown in Fig. 1 still lead to a metalens with a Strehl ratio above 0.9. Such results confirm that the accuracy of our model compared to FDTD simulations is sufficient to represent metasurfaces with performance limited by nanofabrication errors. Therefore, our approach is relevant to computing the properties of the nanostructures and shows the interest of a semi-analytical model to perform a quick initial and optimized design.

Optical design with metasurfaces
To effectively integrate metasurfaces into the design process, it is of utmost importance that the designer can modify and optimize its characteristics just like any other optical component of the system. Therefore, the model's implementation in the design process was carried out with this objective as a priority. In this section, a design methodology is proposed and then applied to present "meta-refractive" hybrid optical systems as a demonstration. Figure 2: Flowchart of the ray-tracing implementation in Zemax using a custom DLL to compute the nanocylinder's distribution, the output wavefront and the incident rays' refraction angle.

Design process
The optical design software Zemax has been selected for our first proof of concept. Its ability to process user-defined surfaces via a dynamic link library (DLL) is an adequate solution to handle calculations related to the semi-analytical model. Its operating principle is shown schematically in Fig. 2. First, Zemax initiates rays and propagates them until the metasurface is reached. Once intercepted, the DLL starts by defining the metasurface from the characteristics of the meta-atoms and the target phase shift profile (subsection 5.1.1). The model then computes the wavefronts so that the refraction angles of the rays are given as output by the DLL (subsection 5.1.2). Exiting the DLL, traditional ray-tracing is resumed across the rest of the optical system. Ultimately, the lens designer can use performance analysis tools to evaluate the system (subsection 5.1.3).

Defining the metasurface
Directly from the Zemax GUI, a lens designer can choose suitable materials and dimensions for the nanostructures. The first step of the design process starts by defining the physical characteristics of the structures to be used: material, dimensions and spacing between structures. For a given height, it is necessary to define an interval of radius (a min , a max ) allowing an adequate phase covering, i.e. |ϕ(a max ) − ϕ(a min )| ≥ 2π. A lens designer can find such information with simple preliminary analyzes using the model corroborated by numerical simulations.
The design phase shift can be expressed in two ways. First, if one desires to produce metalenses specifically designed for focusing, a hyperbolic phase profile is best suited (Eq. 14). Such phase shift leads to a perfect focal spot on the optical axis at a distance f for a plane incident wavefront. However, for a metasurface located between multiple elements, a polynomial phase profile offers much-needed additional degrees of freedom. The phase relation should be defined as Eq. 15.

Computing the output wavefront
Once the metasurface target phase shift is defined, the model is then put to work. Starting from the design phase profile ϕ(r, λ d ), the nanorods' distribution a(r, λ d ) that reproduces exactly ϕ(r) at the design wavelength λ d is found. To do so, the phase relation developed in section 3 should be inverted in order to obtain a distribution of propagation constant as β = ϕ(r)/H. Since the only unknown of this relation is now the radius a of the nanocylinder, the radius distribution needed to obtain a target phase distribution is computed.
For the sake of this numerical application, a continuous distribution of structures of varying dimensions is assumed (see inset of Fig. 2). From the spatial distribution of the structures, the model is again used to determine the resulting wavefront for each wavelength and field of view. Finally, the ray-tracing equations developed previously (Eq. 3 to 5) define the angles of refraction of the incident beam from the phase gradient of the metasurface. Finally, the ray's updated information is given back to Zemax so that ray tracing can continue throughout the system.

Coming up with a design
One of the advantages of our approach is to be able to exploit directly the optimization algorithms included in the design software. Thus, it is possible to optimize for the most efficient optical system at several wavelengths and fields of view, with a proper merit function, while considering the unique behaviour of metasurfaces. Moreover, a lens designer can study how a slight variation in the dimensions of the structures can impact the performance and may enhance the system.
The optimization completed and the performance deemed sufficient, one might want to perform a confirmation with full-wave simulations. The wavefront obtained by simulations, possibly slightly different from the one predicted by the model, can be substituted in Zemax ® . If necessary, fine adjustments can be made to downstream elements to compensate for these differences and regain the target optical performances. This way, optimizing performances using the model in the first place provides an adequately accurate starting point that only requires realistic adjustments when moving to numerical FDTD results.
Ultimately, the resulting metasurface is fabricated using conventional nanofabrication methods. Specialized software such as MetaOptics [8] makes it possible to obtain a technical drawing of the distribution of nanostructures to be designed.

Hybrid refractive and metalens designs
To demonstrate the efficiency and versatility of our method for obtaining metasurface-based optical systems, we now present several variations of a triplet composed of two refractive elements and a metasurface (see Fig. 3). The decision to design a triplet was taken considering that this type of system is the simplest that can correct all of the seven Seidel aberrations over a wide field of view [15].
According to our process, the first step is defining the structures used in the design. For continuity purposes, we have chosen the nanocylinders presented in section 4 for the metasurface, knowing that these respect the phase coverage condition and that the model agrees with numerical FDTD simulations. We then arranged the metasurface as either the last (Fig. 3a), second (Fig. 3b) or first (Fig. 3c) element among a combination of two lenses. The arbitrary goal being to design an imaging system with an f /# = 6 suitable for a ±10 • field of view (FoV), we optimized the metasurface's phase shift profile in addition to the characteristics of the other lenses to obtain a system with almost diffraction-limited performance. For reasons mentioned earlier, we chose a polynomial phase profile over its hyperbolic counterpart to achieve this task. All three systems display good optical performances, a fact suggested by the focal spots of a smaller size than 10 microns. The Airy spot being 8.6 microns in diameter for these systems means that almost diffraction-limited performance is achieved for the whole field of view. An analysis of the modulation transfer function (MTF) supports these conclusions. Indeed, high contrast is achieved at resolutions suitable for most imaging requirements. As mentioned previously, we can expect a slight variation of the predicted phase compared to the FDTD simulations generating a degradation in the optical quality of the metasurface. However, we have determined that for errors listed in Fig. 2, the greatest degradation on the RMS spot size was 0.1 µm for system 3a), 0.2 µm for 3b) and 0.4 µm for 3c). Such variations do not represent a considerable loss of optical quality, highlighting the adequate accuracy of our model. The optimization process made it possible to adapt the metasurface, just like conventional optics, to the system and thus lead to similar overall performance regardless of its positioning. It demonstrates one of the main advantages of our method: considering metasurfaces as a standard optical element in the optimization process while obtaining an accurate performance evaluation.

Conclusion
To conclude, we present a new method to implement metasurfaces in the lens design process via a semi-analytical model used in conjunction with ray-tracing software. We had to adapt the ray tracing equations to consider a phase shift induced by metasurfaces in the first place. Instead of going through lengthy numerical simulations to predict the phase shift, the latter is provided by a semi-analytical model that we have developed. Based on an effective medium theory, the Fabry-Perot effect and a closed-form approximation of the coupled-wave theory, our model predicts the phase shift introduced by cylindrical nanostructures. We then proceeded to its implementation in the optical design software Zemax through a DLL. Starting from a target phase shift for the metasurface, the model integrated into the DLL renders the distribution of nanostructures and calculates the output wavefront so that a ray of a given wavelength and orientation can continue through the system. It led to the development of a 3-step feedback process exploiting the ability to include the metasurface in the optimization routine of an optical system. For demonstration purposes, we have designed hybrid refractive and metasurface systems composed of two refractive elements and a metasurface. The latter, being located at distinct positions for each system, was able to adapt just like conventional elements to render a system with near-diffraction-limited performance across the entire field of view.
Further developments in the generalization of the method to other types of meta-atoms are to be expected. For example, using intrinsically achromatic structures or efficient structures with a large field of view would make it possible to exploit the full potential of metasurfaces in optical design. Ultimately, experimental work could corroborate the accuracy of the semi-analytical model and the relevance of our approach.

APPENDIX A : CMT for nanocylinders
The sum of mode for the radial component inside (1) and outside (2) of the waveguide can be expressed as [27] E (1) r (r, θ) = 1 p 2 1 −iβp 1 A n J n (p 1 r) + ωµ 0 n r B n J n (p 1 r) e inθ (16) −iβq 2 C n K n (q 2 r) + ωµ 0 n r D n K n (q 2 r) e inθ , where A n , B n , C n and D n are arbitrary constants. However, to be solutions of the dispersion relation of the waveguide, the constants should follow these relations where a is the radius of the nanocylinder under consideration. If we only considering the fundamental mode HE 11 , the θ dependency is not present. Thus, the surface integral in one dimension can be computed with the following equations, where [f (x)] * is the complex conjugate of the function f (x).