# Modeling of centimeter-scale metasurfaces in imaging systems

Enzo Isnard, Sébastien Héron, Mahmoud Elsawy, Stéphane Lanteri

► To cite this version:

Enzo Isnard, Sébastien Héron, Mahmoud Elsawy, Stéphane Lanteri. Modeling of centimeter-scale metasurfaces in imaging systems. Computational Optics 2024, SPIE, Apr 2024, Strasbourg, France. pp.15, 10.1117/12.3015928. hal-04660338

HAL Id: hal-04660338

https://hal.science/hal-04660338v1

Submitted on 23 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling of centimeter-scale metasurfaces in imaging systems

Enzo Isnard<sup>a,b</sup>, Sébastien Héron<sup>b</sup>, Mahmoud Elsawy<sup>a</sup>, and Stéphane Lanteri<sup>a</sup>

<sup>a</sup>Université Côte d'Azur, Inria, CNRS, LJAD, 06902 Sophia Antipolis Cedex, France

<sup>b</sup>Thales Research & Technology, 1 avenue Augustin Fresnel, 91120 Palaiseau, France

## ABSTRACT

In recent years, metasurfaces have emerged as promising components for reducing the size and complexity of optical imaging systems. Nevertheless, their incorporation into these systems remains a challenge due to the difference between the numerical tools used in the design of conventional imaging systems (ray tracing) and the design of metasurfaces (full-wave solvers). In this paper, we describe a methodology that combines the local periodicity approximation and a free space propagation technique to rapidly compute the optical response of a large-scale metasurface. We then explain the coupling of this fast solver with a ray tracer to accurately simulate an optical system containing a metasurface.

**Keywords:** Metasurfaces, Imaging systems

## 1. INTRODUCTION

Metasurfaces are 2D optical components structured at sub-wavelength scale that locally control the properties of incident light. They can perform various functions such as beam steering, polarization control and focusing.<sup>1–3</sup> The degrees of freedom they offer and their compactness make them promising components for reducing the size and complexity of imaging systems.<sup>4,5</sup>

However, metasurfaces are still difficult to incorporate into imaging systems due to the difficulty of modeling their behavior together with other components. For a traditional optical system composed of mirrors or refractive lenses, ray-tracing tools are used to predict the imaging performances. For metasurfaces, we can't use these tools, as geometrical optics is no longer valid for describing the interactions of light with sub-wavelength elements. To accurately simulate them, Maxwell's equations have to be solved using finite-differences or finite-elements methods. These methods are computationally expansive and are only adapted for components of around ten wavelengths. Thus, they are of no practical use to simulate metasurfaces inside imaging optical systems as they need to be in the centimeter scale to have industrial applications in imaging. Beyond this size limitation, one may need to put a metasurface neither at the entrance nor at the output of the system. This positioning often leads to smaller incident angle and/or to decrease of the component area because of its relative position with respect to the system stop and pupils. Besides, this also leads to modelling drawbacks: a) the incident fields are not plane waves anymore, and b) computed electromagnetic fields after the metasurface need to be recast as rays.

To solve these issues, we first describe a methodology to compute rapidly the optical response of a meta-surface to an incident field with a non-planar wavefront. We then explain how to couple this fast solver with a ray-tracer in order to simulate a whole system containing mesurfaces and refractive components.
<figure><img src="image_2.png" /><figcaption>Figure 1. Example of an imaging system containing a metasurface. To simulate the whole system, we need to convert rays into electromagnetic (EM) fields and vice versa. As the metasurface lies in the middle of the system, the incident wavefront is not necessarily planar.</figcaption></figure>

## 2. MODELLING OF LARGE SCALE METASURFACES

This part focuses on the methodology to compute the optical response of a single metasurface component. To overcome the difficulty of simulating a large-scale metalens ($10^3 - 10^4 \lambda$), we must approximate the response of metasurfaces with as less computational complexity as possible. We choose to use the local periodicity approximation (LPA) because of its simplicity as it characterizes the optical response of an unit-cell with only one complex transmission coefficient if its size is less than the operating wavelength.⁶ The transmission coefficient can be easily obtained for a continuous range of wavelengths, incident angles and geometrical parameters through the use of cubic spline interpolation for instance. This approximation is sufficiently accurate provided that we use unit cells with low near-field coupling and if the optical power of the system is mainly provided by the refractive elements to limit geometric variations between adjacent unit-cells of the metasurface. In our work we use cylindrical meta-atoms in $TiO_2$ on a $SiO_2$ substrate as an illustration. In Fig. 2 we provide the geometrical parameters of the unit-cells together with the phase-shift and transmission obtained at $\lambda = 600$ nm for different angles of incidence. All simulations were performed using our in-house full-wave solver.⁷

<figure><img src="image_3.png" /><figcaption>Figure 2. In (a) we provide the geometrical shape and its associated parameters of the unit-cell we use. In (b) and (c) we show the phase and transmission response of the meta-atoms for various angle of incidence ($\theta_{inc} \in [0^\circ, 16^\circ]$) at $\lambda = 600$ nm</figcaption></figure>

In order to quantify the validity of the LPA we compare the Point Spread Function (PSF) of a $f/3$ hyperbolic metalens with a diameter of $15 \ \mu m$ illuminated with a $x$ polarized plane wave obtained with a full-wave simulation and with the LPA. First we use the Stratton-Chu formula⁸ to propagate the near-fields above the metasurface to the focal plane. We observe a good agreement between the two PSFs as we have only 3.6% relative error between the two. We also compare full-wave results with LPA when retaining only the main vectorial component from
the incident field and using the scalar Rayleigh-Sommerfeld (RS) formula as near-to-far field transform. The two diffractions formula have only 1.7% relative difference, showing an extremely small variation of Strehl ratio that does not prevent accurate simulation of system performance.

<figure><img src="image_4.png" /><figcaption>Figure 3. Comparison between the PSFs obtained with using the LPA and full-wave results for a f/3 metalens with a diameter of 15 µm.</figcaption></figure>

Instead of using the RS formula we use the Angular Spectrum Propagation method (ASPM) instead, as the two formulas are equivalent⁹ and the second one allows us to use the FFT algorithm to reduce the computational cost of free space propagation:

$$E(x, y, z) = \mathcal{F}^{-1}[\mathcal{F}[E(x, y, 0)]H(\xi, \eta))]$$

where $\mathcal{F}$ is the Fourier transform operator and $H(\xi, \eta) = \exp(i2\pi z\sqrt{\lambda^{-2} - \xi^2 - \eta^2})$ is the transfer function of free space. In the Fourier transform computations, we use zeros padding and apply a low pass filter to avoid aliasing effects.¹⁰ For the computations of Fourier transforms, use of parallel computation drastically reduce the simulations time. In Fig.4, we show a comparison of simulation time between the two and observe that using parallel computing speeds up calculations by up to a factor 6. The limitation of the size of the metasurface we can simulate is due to the memory requirement of saving the sampled near-field in an array. For example, a Graphical Processor Unit (GPU) with 6 GiB of memory can simulate a metasurface of 2.25 mm, which represents 7500² unit-cells.

<figure><img src="image_5.png" /><figcaption>Figure 4. In (a) we show the computations times using a GPU (Nvidia Geforce RTX 2060 GPU) and a CPU to propagate the near-field obtained with the local periodicity approximation using the angular spectrum propagation method as functions of $D$, the diameter of the metalens simulated. In (b) we show the memory consumed by the angular spectrum propagation method on a GPU.</figcaption></figure>
We summarize the accuracy and computational time of each method simulate a $f/3$ hyperbolic metalens with a diameter of 15 µm in Tab.1. The field is calculated on a 201 × 201 grid for each method.

<table><caption>Table 1. Computation time and relative differences with full-wave results of each method. The field is calculated on a 201 × 201 grid for each method using a GPU (Nvidia Geforce RTX 2060 GPU). The diameter of the metalens simulated is 15 µm.</caption><thead><tr><td>Method</td><td>Computational time</td><td>Relative difference with Full-wave</td></tr></thead><tbody><tr><td>LPA Stratton-Chu</td><td>87s</td><td>3.6%</td></tr><tr><td>LPA Scalar RS</td><td>77s</td><td>5.2%</td></tr><tr><td>LPA ASMP</td><td>0.2s</td><td>6%</td></tr></tbody></table>

In conclusion we show that the LPA together with a near-to-far-field transform is an efficient way to obtain the optical response of a large-scale metasurface. We end up using the ASMP method as propagation method as it is by far the most efficient method in term of computational cost.

### 3. COUPLING WITH A RAY TRACER

The hybridization of the proposed methodology with a ray-tracer also requires that we provide an input field to the solver and that we convert resulting fields to rays. For the first step we employ the field tracing equations¹¹ to obtain a field from ray bundles that are traced through the optical system. A ray bundle consists of a main ray and three adjacent rays, which are used to measure the geometric spread as the rays propagate through the system. Given a ray bundle on an input plane associated with a field value $E_{in}$, the value of the field after the ray propagation through the system to an output plane is given by the following formula:

$$E_{out} = \sqrt{\frac{A_{in}}{A_{out}}} e^{i \frac{2\pi}{x} OPL} E_{in}$$

where $A_{in}$ and $A_{out}$ are the area of the triangles formed by adjacent rays on the input plane and output plane respectively and OPL is the optical path length travelled by the main ray as it passes through the system from the input plane to the output plane.

<figure><img src="image_6.png" /><figcaption>Figure 5. Illustration of the field tracing procedure in 2D. The blue line represents a main ray and the green dashed lines its adjacent rays.</figcaption></figure>

Once we computed the incident field with the ray-tracer we can compute the transmitted near-field of the metasurface by approximating locally the incident field to meta-atoms by plane waves. We can then apply the same methodology as in Sec.2 but using the transmission coefficient associated with the incident angle of the ray hitting the center of meta-atoms. With this method we can generalize the LPA to non-planar incident wavefronts that have small curvature compared to meta-atoms size.
To convert the obtained near-field to ray bundles we exploit the fact that rays are colinear to the phase gradient. The phase gradient cannot be obtained directly from the near-field of the metasurface because this field is discontinuous when we approximate it with the LPA. Thus, it is necessary to propagate it over a short distance (typically a few wavelengths) to obtain a smooth wavefront for which we can compute the phase gradient. After this step, we calculate the phase gradient in either the x direction or y direction using the expression of the derivative of the phase with respect to the field:

$$ \partial_x \varphi = \frac{\Re(\partial_x E)\Im(E) - \Im(\partial_x E)\Re(E)}{|E|^2} $$

and we compute $\partial_x E$ using the fact that $\partial_x E = 2\pi i \mathcal{F}^{-1}[\xi \mathcal{F}[E]]$.

The proposed methodology allows us to simulate a metasurface inside an optical system with a low computational cost. The additional rays required to obtain the amplitude of the field incident on the metasurfaces have no significant impact on simulation time, since ray-tracing calculations are significantly less expensive than ASMP. The total computational cost of performing all the calculations on a GPU is 8 seconds.

## 4. CONCLUSION

The methodology we propose enables us to simulate optical systems that combine refractive optical elements with metasurfaces. It consists in calculating the incident field on metasurfaces inside an optical system from ray bundles using the field tracing equations and approximating the local response of the metasurface's constitutive elements using the local periodicity approximation. We generalize local periodicity approximation to non-planar wavefronts that have low curvature relative to the size of the meta-atoms by locally approximating them with plane waves. The transmitted near-field of the metasurfaces is efficiently calculated using the angular spectrum propagation method. The resulting field is then converted into rays using the fact that they are collinear with the phase gradient, which is calculated using the spectral derivative of the fields. With the proposed methodology, we are able to calculate the optical response of a $f/3$ metalens with a diameter of 2.25 mm operating at $\lambda = 600$ nm with near-full-wave accuracy using a GPU with 6GiB of memory. The entire process of calculating the near-field transmitted by the metalens and converting it into rays takes just 8 seconds.

## ACKNOWLEDGMENTS

This work was supported by the French government, through the UCAJEDI Investments in the Future project managed by the National Research Agency (ANR) under reference number ANR-15-IDEX-01. The authors are grateful to the OPAL infrastructure and the Université Côte d'Azur's Center for High-Performance Computing for providing resources and support. We gratefully acknowledge the French National Association for Research and Technology (ANRT, CIFRE grant number 2022/0739) for their financial support.

## REFERENCES

[1] Genevet, P., Capasso, F., Aieta, F., Khorasaninejad, M., and Devlin, R., "Recent advances in planar optics: from plasmonic to dielectric metasurfaces," *Optica* 4(1), 139-152 (2017).

[2] Zaidi, A., Rubin, N. A., Dorrah, A. H., Park, J.-S., and Capasso, F., "Generalized polarization transformations with metasurfaces," *Optics Express* 29(24), 39065-39078 (2021).

[3] Kim, I., Martins, R. J., Jang, J., Badloe, T., Khadir, S., Jung, H.-Y., Kim, H., Kim, J., Genevet, P., and Rho, J., "Nanophotonics for light detection and ranging technology," *Nature nanotechnology* 16(5), 508-524 (2021).

[4] Lee, D., Gwak, J., Badloe, T., Palomba, S., and Rho, J., "Metasurfaces-based imaging and applications: from miniaturized optical components to functional imaging platforms," *Nanoscale Advances* 2(2), 605-625 (2020).

[5] Billuart, J., Héron, S., Loiseaux, B., Amra, C., and Lequime, M., "Towards a metasurface adapted to hyperspectral imaging applications: from subwavelength design to definition of optical properties," *Optics Express* 29(21), 32764-32777 (2021).
[6] Pestourie, R., Pérez-Arancibia, C., Lin, Z., Shin, W., Capasso, F., and Johnson, S. G., “Inverse design of large-area metasurfaces,” *Optics express* **26**(26), 33732–33747 (2018).

[7] “DIOGENeS: A discontinuous-galerkin based software suite for nano-optics.” https://diogenes.inria.fr.

[8] Török, P., Munro, P. R., and Kriezis, E. E., “Rigorous near-to far-field transformation for vectorial diffraction calculations and its numerical implementation,” *JOSA A* **23**(3), 713–722 (2006).

[9] Goodman, J. W., [*Introduction to Fourier optics*], Roberts and Company publishers (2005).

[10] Matsushima, K. and Shimobaba, T., “Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields,” *Optics express* **17**(22), 19662–19673 (2009).

[11] Wyrowski, F. and Kuhn, M., “Introduction to field tracing,” *Journal of Modern Optics* **58**(5-6), 449–466 (2011).