Simple strategy for simulation of large area of axially symmetric metasurfaces

Metalenses are composed of nanostructures for focusing light and have been widely explored in many exciting applications. However, their expanding dimensions pose simulation challenges. We propose a method to simulate metalenses in a timely manner using vectorial wave and ray tracing models. We sample the metalens' radial phase gradient and locally approximate it by a linear phase response. Each sampling point is modeled as a binary blazed grating, employing the chosen nanostructure, to build a transfer function set. The metalens transmission or reflection is then obtained by applying the corresponding transfer function to the incoming field on the regions surrounding each sampling point. Fourier optics is used to calculate the scattered fields under arbitrary illumination for the vectorial wave method and a Monte Carlo algorithm is used in the ray tracing formalism. We validated our method against finite difference time domain simulations at 632 nm and we were able to simulate metalenses larger than 3000lambda0 in diameter on a personal computer.


Introduction
Metalenses consist of subwavelength nanostructures that locally modify the phase profile of an incoming beam to focus it [1][2][3].The unprecedented degree of light manipulation, compactness, and compatibility with standard nanofabrication processes make them attractive substitutes for conventional refractive optics systems.These cutting-edge devices have been designed for a wide range of wavelengths (from 50 nm to 3600 nm) and are now developed for industrial applications [3].Diffraction limited focusing [4], wide field of view imaging [5,6], achromatic focusing [7], and endoscopic imaging [8] are a few examples of the breadth of applications enabled by metalenses.
Recently, mass-manufacturable metalenses with diameters on the order of a few centimeters have been demonstrated [9][10][11][12].Normally, metalenses are simulated under certain approximations or advanced techniques that reduce the computational burden.Conventional rigorous numerical methods to solve Maxwell's equations, such as finite difference time domain (FDTD) and finite elements method (FEM), are not suitable because they require enormous computational resources [13].Several strategies to accelerate these methods have been proposed, including hardware acceleration for the FDTD method [14], augmented partial factorization (APF) [15], and low-overhead distribution on a GPU-based simulation [16].These approaches have successfully simulated metalenses with dimensions as large as 600 0 [16].However, the most common approach to simulate even larger area metalenses is the local approximation method [6,14,17,18], where the transmitted field by each nanopost is assumed to be constant and equal to its array response.This approach, however, does not account for the coupling among nanostructures [1] and cannot simulate metagratings .
Here, we propose a method to accurately simulate large area metalenses in a timely manner.Our strategy is based on sampling the metalens phase gradient profile and modelling each point as a binary blazed grating using the nanopost or metagrating design of choice.We then build a transfer function library for each blazed grating under plane wave incidence for different incident angles and polarizations.A similar approach has been recently demonstrated to simulate quantum emitters' response near periodic patterned hyperbolic metamaterials [19,20].Our model allows us to simulate the metalens response to an arbitrary field distribution under a wave vector and ray tracing models.The model expands the field and rays in terms of the diffraction orders of the metalens, allowing for an unprecedented analysis of the focused field.We compared our method against rigorous FDTD simulations and managed to reduce the simulation time and memory requirements by at least on order of magnitude whilst keeping the resulst accurate.We used our method to simulate metalenses with diameters larger than 3000 0 and separate the field into the contributions of each diffraction order.

Metalens model rationale
We propose a model based on a transfer-function approach that considers coupling among different posts.The metasurface phase gradient is sampled on  positions, as illustrated in Fig. 1.The   ℎ patch can be modelled as a blazed binary grating with period given by the grating equation as where  (ì  1 ) ≡  1 = ()/ is the radial phase gradient value calculated at the   ℎ sampling point.Any radial phase profile can be sampled using our phase gradient library, as shown in Fig. 1.The blazed binary grating is modeled according to the nanostructure or metagrating design being used such as rectangular, elliptical posts, among others, as shown in Fig. 1 (c).This approximation is solely used to obtain the transmission and/or reflection coefficients of that region and do not account for the phase profile curvature.We address this issue by locally correcting the phase gradient according to the ideal phase profile, as will be discussed later.We calculate the blazed binary grating transmission and reflection coefficients under plane wave excitation at different angles of incidence and polarizations using the rigorous coupled-wave analysis (RCWA) method [21,22].The transfer function is defined as the diffraction efficiencies or field amplitudes of the transmitted/reflected waves as a function of the incoming in-plane wave-vector.We solve Maxwell's equations for a supercell of the blazed binary grating on a region with a cross-section area given by  ×   , where  is the unit cell size of a single post on the metalens design, as shown in Fig. 1 (c).Note that depending on the metasurface nanopost unit cell size, we may have to use a supercell larger than   to fit an integer number of posts inside it.With this approach, the coupling among the posts is fully modeled, but the coupling between different regions is neglected.To generate each blazed binary grating, we perform the phase profile sampling as necessary.After calculating the transfer functions, they are used to obtain the scattering properties of the metalens.The transfer function calculation needs to be realized once, and it can be hastened using accurate and faster methods [14][15][16].To save computation effort, we calculate a single transfer function for each radial sector to later rotate them to the new system of coordinates taking advantage of the metalens rotation symmetry.

Vectorial wave model
The metalens is modeled as a spatially patched transfer function that modulates an arbitrary field distribution and relies on the angular spectrum formalism for the propagation in free space [23].As explained in the previous section, we sample the metalens phase profile gradient with linear patches as shown in Fig. 2(a).Moreover, we also split the metalens along the azimuthal direction, leveraging its rotation symmetry to use the same transfer function but properly spatially rotated, as shown in Figs.2(b) and 2(c).That is, we approximate the metalens phase profile by a piece-wise linear phase profile at the points ì  ,  creating  sectors radially, and we split each sector into   regions where  indicates a given radial sector, as represented in Fig. 2. The total transmitted or reflected fields by the metalens are calculated as the sum of the contribution of each sector where ì  ,  (ì ) is the electric field distribution calculated on sector  of region  and    (ì ) is a window function, as shown in Fig. 2(c), that limits the region area and is given by Eq. (S2).ì  is defined as a radial vector centered on the metalens (see Fig. 2 (a)) .If the incoming shows the angular splitting for each sector, forming a given region.ì  is defined as the position vector centered on the metalens and ì    is the vector pointing to the center of the region ĳ.ρ  and φ  are the blazed direction versor and the corresponding orthogonal versor of region ĳ, respectively.field distribution has a Fourier transform given by ì  0 ( ì  ∥ ) (bold symbols represent the Fourier transformed quantities), the transmitted field spatial spectrum on axial region j of sector i (see Fig. 2) is given by where is the transfer function defined as a tensor given by Eq. (S4).Eq. ( 3) assumes that the metalens acts as a linear operator on the incoming field.Such ansatz is corroborated by the linear property of Maxwell's equations [19,20,23].As shown in the SI, after operating the transfer function on the incoming field distribution, the transmitted field spectrum can be calculated as where ì    ( ì  ∥ ) is the spectrum of the output field produced by the -th diffraction order and is given by Eq. (S10).The reciprocal vectors are defined according to the versors ρ  and φ  , shown in Fig. 2 (c).Finally, to obtain the field in real space, we simply perform an inverse Fourier transform in ì    ( ì  ∥ ).Applying the inverse Fourier transform and its shifting property in Eq. ( 4) we have that As seen in Eq. 5, the total field on the (, ) region is given by the linear superposition of the diffracted fields modulated by the linear phase of the order.Furthermore, the phase term ∥ controls the central position of the diffracted field spectrum and creates a local linear phase profile distribution, which is a consequence from the linear patching approach.This effect induces aberrations on the wavefront and can be detrimental to the properties of the metalens.To eliminate this problem, we correct this phase term by substituting it to the original phase profile, (ì  ∥ ), as discussed in the SI.

Ray tracing
Although more precise and faster in describing the behavior of large-area metasurfaces, the computational burden of the vectorial approach can be further reduced by applying the proposed method on the ray tracing method.The metalens is now treated as a phase discontinuity on the ray path [24] but also accounting for the diffraction efficiency of the metalens.Given a ray with in-plane momentum ì  ∥ incoming at a point ì  ∥ on the metalens (see Fig. 2(a)), we can model the diffraction by treating it as a probabilistic event with a Monte-Carlo algorithm.That is, there are three possible outcomes for the ray in our model after it interacts with the metalens: it can either be diffracted back (reflection), forward (transmission), or absorbed if the structure is lossy.The probabilities are taken as the diffraction efficiencies of the patch approximation obtained by solving Maxwell's equations.For the scattering processes we define   ( ì  ∥ , ì ) and   ( ì  ∥ , ì ) as being the probabilities of a ray with incoming in-plane momentum ì  ∥ and polarization state ì  being diffracted in transmission and reflection, respectively, and  ∈ D is a given diffraction order in the set of stored diffraction orders space (D).From energy conservation, we can define the probability of a ray being absorbed as , as shown in the SI.
We define a cumulative probability distribution  as where, , ] and we omitted the region indices for clarity.Therefore, we can use  to define intervals where a given outcome might happen.The diffraction order can then be found by generating a random number and analysing in which interval it has fallen.Given the uniformly distributed random variable  ∈ [0, 1], we can obtain , and consequently the outcome of the event from the sequence , by solving the following equation where  is the resulting diffraction index and index of the sequence , which is used to map a given event.With the scattering event determined, the resulting ray can either be scattered or absorbed.If it is absorbed, then no other ray is produced.On the contrary, if the ray is reflected or transmitted, then a new ray is generated with momentum given by where n is the normal vector to the metasurface, ì   is the reciprocal vector of the g-th order, and   is the resulting ray wave-vector along the normal direction, which can be found from the dispersion equation in the medium.Note that the reciprocal vector radial component is corrected locally according to the local phase gradient instead of using the blazed binary grating momentum.As shown in the SI, we can also obtain the diffracted polarization by the metalens.Each ray can be labelled according to the corresponding diffraction order it originated from, allowing a better understanding of the metalens behavior, without the need to perform separate simulations, see Figs. (S1) and (S2) in the SI for more detail.Such discrimination is also possible on the vectorial wave method, but it would require storing a field distribution of each order, which is impractical.We compare the focusing profile of a metalens operating at 632 nm against FDTD simulations.The ray tracing map distribution is obtained by calculating the ray density crossing the plane  = 0 and the field distributions were normalized with respect to their peak values on the focal plane to highlight the field distribution.The simulation parameters are discussed in section 3.3.We use 1200 nm tall glass-based nanoposts in our design [12,25].The metalens focal length is 200 m with a diameter of 100 m (NA = 0.25) and it encodes a hyperbolic phase profile.The vectorial wave based model qualitatively reproduces the FDTD field distribution with good fidelity, even accounting for the appearance of higher-order focal spots as shown in Figs.3(a) and 3(b), respectively, at normal incidence.The model can also reproduce the field distribution at 20  of incidence as shown by Figs.3(e) and 3(f), respectively.The ray tracing model also simulates the focusing field distribution and the high order focal spots as shown in Figs.3(c) and 3(g).However, it accounts only for the power flow of the diffracted rays and idealizes the intensity distribution as it does not model interference among the rays.The interference leads to the well-known diffraction limit in optics and can be visualized by the field intensity distribution on the focal plane, as shown in Figs.3(d) and 3(h) at normal and oblique ( 20 • )incidence, respectively.All plots were normalized with to their own power flux and peak intensity of the ideal Airy disk distributrion.The field distributions at the focal plane obtained by our method presents good agreement with the corresponding distributions FDTD at normal and oblique incidences.The transversal cuts of the point spread function at oblique incidence are highly distorted due to off-axis aberration -mainly coma, as shown in Fig. 3(j).We also calculated the focusing efficiency on an circle with diameter ten times larger then the PSF full width at half maximum (FWHM).At normal incidence, we obtained 73.5%, 72.3% and 71.5% with the FDTD, vectorial wave and ray tracing methods respectively.

Focusing profile comparison
Finally, the ray tracing model allows us to easily tag each ray according to the diffraction order it originated.Fig. (4)(a) shows the ray tracing of the simulated metalens at normal incidence.To avoid overcrowding, the plot is limited to 250 rays in total and discriminates the diffraction orders according to the ray color, which are limited to the [−5, 5] range.The negative orders gives rise to additional short real focal spots, whereas the positive ones to virtual focal spots.We can calculate the diffraction efficiencies of each order by tallying the number of rays directed into them.Fig. 4(b) shows the resulting normalized diffraction orders histograms produced at normal and oblique incidence, respectively.We also applied our model to simulate the performance of a doublet metalens based on [26], as shown in Fig. S1 of the SI.

Diffraction efficiency calculation
After qualitatively showing that our models can calculate the focusing field profile of a metalens, we assess the efficiency prediction of our model against numerical simulations performed using FDTD.We used a glass-based nanopost design described in [25] where metagratings are used to increase the metalens numerical aperture.Here, we scan a beam across a metalens and calculate the diffraction efficiency of the first few orders as a function of its position from the metalens center.The metalens focal length is 7 mm with a diameter of 15 mm.We use an unpolarized Gaussian beam with 100 m (158 ) waist.The simulation region for each scan point is 200 m × 200 m.
In particular, the efficiency of the  = −1 order is the highest at normal incidence, and it is almost constant throughout the whole scan, remaining higher than 60%, which is an good indication that the proposed design focuses light efficiently.These results highlight how accurately our model can simulate the metalens even at oblique incidence when coupling among adjacent posts takes over.These results have also been compared against experimental data with a good match, see [25].

Computational resources
In this section, we compare each method's simulation time and memory requirements against FDTD simulations.Note that this estimation does not account for the computational resources used in calculating the transfer functions.We simulate fused silica glass-based metalenses, operating at 632 nm [25], with NA = 0.5 and different focal lengths (f).The FDTD simulations are performed on the commercial software Lumerical inc.using a uniform mesh with 20 nm sampling in all directions.The simulation volume is set to 2R × 2R × 2 , where R is the metalens radius ( =  tan(asin( )).Moreover, the FDTD simulations are performed only to calculate the transmitted near field by the metalens, and the free space propagations can be obtained using the angular spectrum formalism [23].The FDTD simulations are carried out on the FAS cluster with 128 cpus distributed over 4 nodes.The proposed approaches are executed in serial on a personal laptop and have room for improvement if a parallelization scheme is used.One could easily distribute the simulation of each patch to different CPUs to hasten the simulation.The field matrices on the wave optics model are calculated using 200 nm of sampling.Finally, the ray tracing model is calculated using a ray density of 5000 rays/ 2 hitting the metalens.Figs.6(a) and 6(b) show the memory required and the time consumed for each method as a function of the metalens radius (increasing focal lengths).Note that even though the models developed here run in serial, the time required for the simulation is considerably lower compared to the FDTD simulation.Furthermore, the proposed approach requires at least one order of magnitude less memory.

Large area metalens focusing profile simulations
Our approach's low computational time and minimal memory usage enable the simulation of large-area metasurfaces.Here, we simulate the focusing profile of hyperbolic and quadratic metalenses [3,6] operating at 632 nm, utilizing the same post-design described in the preceding sections.The metalenses focal length are 1 mm with a diameter of 2.26 mm (NA = 0.75), which amounts to ∼ 3755 0 .The metalenses are excited by plane waves at normal and 20 • of incidence.Figs.7(a)-7(d) and 7(e)-7(h) show the longitudinal field intensity distribution focused by the metalenses calculated using the wave optics and ray tracing models, respectively.A metalens with this size would take over 100 hours and require over a petabyte of memory to simulate on a FDTD cluster with over 400 cpus running in parallel, as shown in Fig. 6(b), only to obtain the near field.Our model requires approximately two orders of magnitude less time and memory to simulate the same metalens.
Both models display similar focused field distributions, but the ray tracing model only accounts for the power flow without any interference effects.At normal incidence, the mirror-symmetric focusing profile of the hyperbolic profile is obtained in both models, and both feature a weak second-order focusing around z=500 , as shown in Figs.7(a) and 7(e).When illuminated at 20  of incidence, both models show a strong aberrated focal profile due to coma, and a slightly stronger second-order focusing as shown in Figs.7(b) and 7(f).The focusing field distributions of the quadratic profile at normal incidence, obtained with the vectorial wave optics and ray tracing methods, are shown in Figs.7(c) and 7(g), respectively, and do not have mirror symmetry around the focal plane, which is a manifestation of spherical aberration [6] .The quadratic field distribution remains almost the same at oblique incidence, as shown in Figs.7(d) and 7(h), due to its wider field of view [27].

Conclusion
We propose a strategy to simulate large area metalenses by patching it into smaller parts that can be simulated faster using Maxwell's equations.The patching process is similar to the phase sampling used to design a metasurface profile.However, we sample the phase gradient instead, and a corresponding blazed binary grating models each sampling point.Maxwell's equations are then rigorously solved for each linear piece using a single supercell of the blazed binary grating to obtain a transfer function that is a function of the angle of incidence and polarization.The same transfer function can be reused to model different metalenses or metasurfaces with radial phase profiles for a given metasurface post design.Thus, after the transfer function is obtained, the metalens transmitted or reflected fields can be simulated either using a ray tracing approach or a vectorial wave model, depending on the desired application.We found very good agreement between our model, FDTD simulations, and experiment.This approach reduces the time and memory requirements by orders of magnitude, allowing the simulation of metalenses with diameters larger than 3755 0 on a personal computer.
Funding.Content in the funding section will be generated entirely from details submitted to Prism.Authors may add placeholder text in the manuscript to assess length, but any text added to this section in the manuscript will be replaced during production and will display official funder names along with any grant numbers provided.If additional details about a funder are required, they may be added to the Acknowledgments, even if this duplicates information in the funding section.See the example below in Acknowledgements.For preprint submissions, please include funder names and grant numbers in the manuscript.

Fig. 1 .
Fig. 1.(a) and (b) show the sampling of the gradient of an arbitrary radial phase profile and the region they cover on the phase profile, respectively.(c) depicts the equivalent blazed binary grating approximation on the metalens.

Fig. 2 .
Fig. 2. (a) show how the different sectors are set up based on the gradient sampling.(b) shows the angular splitting for each sector, forming a given region.ì is defined as the position vector centered on the metalens and ì    is the vector pointing to the center of the region ĳ.ρ  and φ  are the blazed direction versor and the corresponding orthogonal versor of region ĳ, respectively.

Fig. 3 .
Fig. 3. Field intensity distribution focused by a metalens with a focal length of 200 m and NA=0.25 calculated using different methods.The operating wavelength is 632 nm in all cases, (a)-(d) [(f)-(i)] show the longitudinal field amplitude distributions focused by the metalens at normal [20 • ] incidence calculated using the vectorial wave linear patch, FDTD, ray tracing with linear patching.(e) and (j) show the transversal cuts of the field distributions on the focal plane at normal and 20 • incidence, respectively.

Fig. 4 .
Fig. 4. Ray tracing of a hyperbolic metalens with  = 200 and NA=0.25.(a) and (c) show the resulting ray tracing for the metalens illuminated at normal and oblique (20  ) incidence, respectively.We only show 250 rays on this plot and limited the diffraction order to ±5.(b) and (d) show the diffraction orders normalized histograms of the transmitted rays for normal and oblique incidences, respectively.

Fig. 5 .
Fig. 5. Diffraction efficiencies of different orders scattered by the metalens when scanned radially by an unpolarized Gaussian beam with 100 m of waist and different angles of incidence along the scanning line.From left to right, the columns show the diffraction efficiency of the -2,-1,0,1, and 2 orders, respectively.The first, second and third rows show, respectively, the results at normal incidence, 10 • and 20 • of incidence.The green, red, and blue lines show the ray tracing model, our vectorial wave model, and the FDTD model results.The operating wavelength is 632 nm.The metalens focal length is 7 mm with a diameter of 15 mm.

Fig. 6 .
Fig. 6.(a) and (b) show, respectively, the memory and time required to simulate glass-based metalenses with different radius () (dots) and the estimated values using polynomial interpolation from the computed data (dashed lines).The FDTD simulations are used only to calculate the transmitted near field and are performed using Lumerical on the FAS cluster with 128 cpus distributed over 4 nodes using parallelization.The other models were run in serial and have room for improvement.All simulations are performed on a normal incident plane wave operating at 632 nm.

Fig. 7 .
Fig. 7. (a) -(d) and (e) -(h) show the longitudinal field intensity distribution focused by the metalenses that were calculated using the wave optics and ray tracing models, respectively.(a) and (b) [(e) and (f)] show the focusing profile of a hyperbolic metalens at normal and 20 • of incidence, respectively.(c) and (d) [(g) and (h)] show the focusing profile of a hyperbolic metalens at normal and 20 • of incidence, respectively.The metalenses focal length are 1 mm with a diameter of 2.26 mm (NA = 0.75) and the operating wavelength is 632 nm.