# Towards a metasurface adapted to hyperspectral imaging applications: from subwavelength design to definition of optical properties

Jules Billuart, Sébastien Héron, Brigitte Loiseaux, Claude Amra, Michel Lequime

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

## To cite this version:

Jules Billuart, Sébastien Héron, Brigitte Loiseaux, Claude Amra, Michel Lequime. Towards a metasurface adapted to hyperspectral imaging applications: from subwavelength design to definition of optical properties. Optics Express, 2021, 29 (21), pp.32764. 10.1364/OE.432969. hal-03454629

HAL Id: hal-03454629

https://hal.science/hal-03454629

Submitted on 29 Nov 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Towards a metasurface adapted to hyperspectral imaging applications: from subwavelength design to definition of optical properties

JULES BILLUART,<fnref n="1" />,* SÉBASTIEN HÉRON,<fnref n="1" /> BRIGITTE LOISEAUX,<fnref n="1" /> CLAUDE AMRA,<fnref n="2" /> ID AND MICHEL LEQUIME<fnref n="2" /> ID

**Abstract:** We numerically demonstrate the capability of a single metasurface to simultaneously separate and focus spectral features in accordance with the specifications of a pushbroom hyperspectral imager. This is achieved through the dispersion engineering of a library of two-level TiO₂ nano-elements. Sommerfeld integrals are used to confirm our numerical simulations provided by our solver based on Fourier modal method. As a proof of concept, a metasurface with a 175 µm diameter is designed to be compatible with hyperspectral imaging over a spectral range of ± 50 nm around 650 nm with a spectral resolution of 8.5 nm and a field of view of 8° around the normal incidence (angular resolution of 0.2°).

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Spectral imaging is an efficient method to determine the electromagnetic properties of objects within a complex scene by capturing light emitted by these objects across several spectral bands for each pixel of the image [1]. Light spectral decomposition can be discrete (multispectral imaging) or continuous (hyperspectral imaging), while the image formation can be achieved pixel by pixel (whiskbroom scanner), line by line (pushbroom scanner) or in a global way (snapshot acquisition). This spectral information is often related to the chemical composition of the imaged objects, which is of great interest in many fields such as food safety [2,3], disease diagnosis [4], or remote sensing [5]. Optical performances of these imaging instruments become more and more demanding, so that the complexity of their design increases, and consequently their size and their mass, which is obviously detrimental for embedded applications. For some applications, this size reduction objective has led to integrate dedicated interference filters (stripe filters, pixelated filters, or linear variable filters [6,7]) as close as possible of the focal plane to simplify the global design of the instrument.

Recently, optical metasurfaces [8], i.e. 2D components spatially structured at the wavelength scale, have emerged as a powerful technological avenue allowing to develop ultra-thin and high quality multifunctional devices, such as polarimetric analyzer [9], optical vortex converter [10], high numerical aperture metalenses [11–14], etc. The purpose of our work is to propose a new design approach to achieve important size gain in the field of hyperspectral imaging.

Performance of a spectral imager is mainly defined by its spatial and spectral resolution, but also by its field of view and spectral range. For classical systems using bulk refractive, reflective, or dispersive components, these features are quantified using ray tracing algorithms and scalar diffraction theory [15]. As the minimum feature size of elements involved in the design of a metasurface becomes closer to the incident wavelength, these standard approaches can no more be used [16] and a full-wave simulation is required to correctly estimate the spatial and spectral dependence of the electromagnetic field in the vicinity of the focal plane [17]. As a matter of fact, Pommet et al. [18] have investigated the limit of use of scalar diffraction theory in the case of

<footnote n="1">¹Thales Research and Technology, France</footnote>

<footnote n="2">²Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France</footnote>

<footnote>*jules.billuart@thalesgroup.com</footnote>
perfectly conducting gratings, and non-negligible errors arise when features size is less than 14 λ [16]. However, such a full-wave simulation is performed at high computational cost, making it inappropriate to simulate macroscopic metasurfaces involving millions of scatterers [19].

Besides, the microscopic features of a metasurface enable the use of new degrees of freedom for the design of a compact optical component. We can play with intrinsic dispersion of the constitutive nano-elements to tailor the dispersive optical properties of the component: this is the purpose of the so-called dispersion engineering [20–25]. This approach, first widely used for the design of achromatic meta-lenses [24,26,27], has recently been extended to the multispectral imaging domain [28,29]. In 2019, Zhu et al. [28] have demonstrated the ability of a single metasurface to separate and focus four discrete wavelengths from the visible spectrum in a focal plane perpendicular to the component. However, the paper does not describe the method used to optimize the design of the proposed configuration, particularly the choice of detector orientation. In the same order of ideas, Pestourie et al. [30] have proposed in 2018 a design allowing to focus three colinear incident beams of different wavelengths in three different points of a plane parallel to the metasurface. The choice of this design has been made by solving an inverse problem [31], and then optimizing the structure with the help of a conservative convex separable approximation variant of the original method-of-moving-asymptotes algorithm (CCSA-MMA) [32].

In this work, we present a method allowing to define and optimize the design of a metasurface able to synthetize all the optical functionalities of a spatial-scanning pushbroom hyperspectral imaging system. This implies that this planar component performs both the dispersion of a continuum of wavelengths over a given spectral bandwidth Δλ and the focusing of all these wavelengths along a focal line, and that for all the incidence angles θ defined by the useful height of the entrance slit of the instrument. In other words, the properties of this metasurface can be described by a functional F taking both input parameters: incident angle θ and wavelength λ. Therefore, Section 2 is devoted to the definition of this functional F, while Section 3 presents the dispersion engineering of the nano-elements involved in the design of this metasurface. After selecting manufacturable configuration of nano-elements, the macroscopic properties of the instrument are then analyzed in Section 4 in terms of radiometry and point spread function (PSF). Its optical performances are then compared to those of an ideal phase mask and a zone plate. Finally, Section 5 illustrates our design methodology on an example characterized by a spectral bandwidth of 100 nm around 650 nm and an angular field of view of 8°, while Section 6 summarizes our main results and provides conclusive remarks.

## 2. Functional adapted to hyperspectral imaging applications

For a pushbroom hyperspectral imaging system, the metasurface, located in the plane z = 0, must ensure a perfect stigmatism between an object point rejected to infinity in a direction θ located in the yOz plane and an image point $F_{θ,λ}$ located in a uv focal plane where v is parallel to y and where the u-coordinate of this image point is a linear function of the wavelength λ (see Fig. 1). This means that the functional $F(x, y, λ, θ)$ that describes the effect of the metasurface must satisfy the following condition:

$$-y \sin \theta + F(x, y, \lambda, \theta) + [MF_{\theta,\lambda}] = C_{\theta,\lambda} \quad (1)$$

where M is the point of the metasurface that is defined by the coordinates x and y, $[MF_{θ,λ}]$ the optical path between M and $F_{θ,λ}$, and $C_{θ,λ}$ a numerical factor that only depends of the angle θ and the wavelength λ. The coordinates of $F_{θ,λ}$ in the u, v, w frame are :

$$[u(\lambda), v(\theta), 0] \quad \text{where} \quad \begin{cases} u(\lambda) = a\lambda + b \\ v(\theta) = c \tan \theta \end{cases} \quad (2)$$
while the set of linear relations allowing to describe the change of system of coordinates is as follows

$$ \begin{cases} x = f_x + u \cos \alpha \\ y = v \\ z = f_z + u \sin \alpha. \end{cases} \qquad (3) $$

Consequently, the functional $\mathcal{F}$ of the metasurface is defined by

$$ \mathcal{F}(x, y, \theta, \lambda) = C_{\theta, \lambda} + y \sin \theta - \sqrt{[x - f_x - (a\lambda + b) \cos \alpha]^2 + [y - c \tan \theta]^2 + [f_z + (a\lambda + b) \sin \alpha]^2} \quad (4) $$

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

<figure><img src="image_4.png" /><figcaption>Fig. 1. Schematic of an optical architecture for hyperspectral imaging system - a) an incident plane wave coming from the left is focused onto a matrix of photodiodes by a metasurface (MS) characterized by a functional $\mathcal{F}$ and an entrance pupil $D_x$ in the $(x, z)$ plane; the detector is tilted by an angle $\alpha$ in this plane - b) sketch of the angular illumination $\theta$ of the metasurface in $(y, z)$ plane; the metasurface has an entrance pupil $D_y$ in this plane. The schematic assumes an entrance slit to keep spatial information along one direction only (y axis). The scan along $x$ axis is provided by the translation of the whole hyperspectral imager along this $x$ axis.</figcaption></figure>

The constant $C_{\theta, \lambda}$ is independent of $x$ and $y$, which implies that

$$ \mathcal{F}(0, 0, \theta, \lambda) = C_{\theta, \lambda} - \sqrt{[f_x + (a\lambda + b) \cos \alpha]^2 + [c \tan \theta]^2 + [f_z + (a\lambda + b) \sin \alpha]^2} \quad (5) $$

By combining the two last equations, we obtain finally

$$ \begin{aligned} \mathcal{F}(x, y, \theta, \lambda) = \mathcal{F}(0, 0, \theta, \lambda) &+ y \sin \theta \\ & - \sqrt{[x - f_x - (a\lambda + b) \cos \alpha]^2 + [y - c \tan \theta]^2 + [f_z + (a\lambda + b) \sin \alpha]^2} \\ & - \sqrt{[f_x + (a\lambda + b) \cos \alpha]^2 + [c \tan \theta]^2 + [f_z + (a\lambda + b) \sin \alpha]^2}. \end{aligned} \quad (6) $$

For the sake of simplicity, we will restrict the optimization of the metasurface design to a unique angle of incidence ($\theta = 0^\circ$), then we will analyze the performance evolution of this spectrally optimized design with respect to $\theta$ in the yOz plane (see Fig. 1(a)). In this restricted case, the phase change induced by the crossing of the metasurface at zero angle of incidence is defined by :

$$ \Phi_{\mathcal{F}}(x, y, \lambda) = \Phi_{\mathcal{F}}(0, 0, \lambda) - \frac{2\pi}{\lambda} \left\{ \sqrt{[x - f_x - (a\lambda + b) \cos \alpha]^2 + y^2 + [f_z + (a\lambda + b) \sin \alpha]^2} - \sqrt{[f_x + (a\lambda + b) \cos \alpha]^2 + [f_z + (a\lambda + b) \sin \alpha]^2} \right\}. \quad (7) $$
### 3. Dispersion engineering

One way to design a metasurface achieving the optical function described in Eq. (7) is to use the direct design approach [33], where the metasurface is considered as an array of independent scatterers, i.e. not interacting with neighbors. The phase profile $\Phi_{\mathcal{F}}(x, y, \lambda)$ is sampled with a period $\Lambda$ smaller than the wavelength and the local phase value is obtained by selecting a nano-element whose effect on the incident wave is in agreement with this value. Moreover, the choice of a subwavelength period ensures the absence of parasitic diffractive orders and provides a very good approximation of the theoretical phase profile.

To illustrate this approach in a simple way in the case of a hyperspectral imaging functionality, the problem is limited to the design of a 1D cylindrical metasurface. For the sake of simplicity, we call a 1D cylindrical metasurface a subwavelength engineered component whose structure and profile are invariant along the y axis (cf. Fig. 2(a)). A 1D problem is chosen to reduce the computational cost while calculating the electromagnetic fields through a full wave simulation; it is important to stress here that one can easily extend this design approach to a 2D configuration. In this 1D version of the design problem, the electromagnetic fields are translation invariant along the y axis, and the spatial phase gradient necessary to achieve the focusing of the incident wave only exists along the x axis (see Fig. 1). For such a cylindrical metalens, the stigmatism condition is only possible along the x axis, because the component does not exhibit any phase gradient along the y axis. Thus, a single point source is focused along the v axis over a length $D_y$ which corresponds to the length of the entrance pupil in the yOz plane (see Fig. 1(a) and (b)). In Section 5.2 is explained how such a cylindrical component is nevertheless able to perform the imaging function associated with different angles of incidence $\theta$.

The 1D nano-elements used in this direct design approach are two-level TiO$_2$ nano-ridges characterized by heights $h_i$ and widths $w_i$ ($i = 1, 2$), as schematized in Fig. 2(a), and located at the top of a fused silica substrate. The choice of titanium dioxide provides a high refractive index contrast between the nano-elements and the surrounding medium (air), which limits the interactions between closest neighbors [11], also known as cross-talk. The phase response of an individual stack of nano-ridges is calculated with the help of Fourier Modal Method (FMM) [34], a version of Rigorous Coupled Wave Analysis [35] implemented in a home-made solver. Each nano-element located at the position $x$ on the metasurface must delay the locally incoming wavefront with an appropriate phase term $\Phi(x, \lambda)$. If the spectral bandwidth $\Delta\lambda$ of the imager contains $N$ elementary wavelengths, each stack of nano-ridge must then fulfill $N$ phase requirements for a given position $x$ along the metasurface. This approach makes the design of the metasurface, even 1D, really challenging. Therefore, it is more convenient to approximate the local phase response $\phi(\lambda)$ by a Taylor expansion up to the first order, i.e.

$$ \phi(\lambda) = \phi(\lambda_d) + (\lambda - \lambda_d) \left. \frac{\partial \phi}{\partial \lambda} \right|_{\lambda_d} + O[(\lambda - \lambda_d)^2] \quad (8) $$

The Taylor expansion of the local phase response $\phi(\lambda)$ has been extensively used to design achromatic metalenses [20,22,26] or any dispersion-engineered metasurfaces [36]. The Taylor expansion up to the second order is commonly used for two reasons. Firstly, it ensures a correct approximation of the nano-element phase response, especially if the bandwidth $\Delta\lambda$ is large. Secondly, for achromatic lenses the spectral phase has a hyperbolic shape, which is even more pronounced as the lens diameter increases [22]. In this case, a simple tangent (namely the first order term) centered at a design wavelength $\lambda_d$ is not sufficient to correctly approximate the phase functional $\Phi_{\mathcal{F}}(x, \lambda)$. In our work, the spectral range $\Delta\lambda = 100$ nm which is small enough to correctly estimate the nano-element spectral response (see Fig. 2(b)). In this work, we only consider the TE polarization (Fig. 2(a)), but the design approach can be obviously applied to a TM polarized field, for which the phase values are obviously no longer the same due to the asymmetric nature of 1D nano-elements. The analysis has been conducted in TM mode to see if
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Dispersion engineering in TE polarization. a) Two-level nano-ridge 1D structure
- b) red dots, spectral dependence of the phase response of a nano-element; blue curve,
linear fit of this spectral dependence at $\lambda_d$ = 650 nm - c) mapping of the phase space
[$\phi(\lambda_d)$, $(\partial\phi/\partial\lambda)_{\lambda_d}$] corresponding to the nano-element library made of single level ridge;
the color code is defined by the right bar and corresponds to the mean transmission $|t|$ over
the design bandwidth - d) mapping of the phase space corresponding to the nano-element
library made of two-level ridge system, showing a denser filling of the phase space - e)
comparative plots of the phase requirements of the metasurface (green dots) and the phase
space covered by the nano-element library (same color code as in graph c).</figcaption></figure>
substantial changes would appear, and the results reported here after remains globally the same as in the TE polarization.

With such a first order Taylor polynomial, the phase requirements for the nano-elements reduce to 2 terms, the phase at the design wavelength $\phi(\lambda_d)$ and the derivative of the phase at the same wavelength $(\partial\phi/\partial\lambda)_{\lambda_d}$. In Fig. 2(b), the phase response is studied in the wavelengths interval from $\lambda_{min} = 600$ nm to $\lambda_{max} = 700$ nm, and the partial derivative is computed at a design wavelength $\lambda_d = 650$ nm corresponding to the middle of this interval. Because of the two-level structure of the nano-ridges, the number C of possible permutations is a function of geometrical characteristics of the unit-cell (see Fig. 2(a)). As an example, for a subwavelength period $\Lambda = 400$ nm, a minimum thickness $e_{min} = 60$ nm, a spacing $s = 30$ nm at the border of each unit-cell, and a manufacturing step $\delta_e = 1$ nm, the possible widths for each level are given by

$$ \begin{cases} w_1 = e_{\min} + n_1 \delta_e, & n_1 = 0, 1, \dots, (\Lambda - s - e_{\min})/\delta_e \\ w_2 = e_{\min} + n_2 \delta_e, & n_2 = 0, 1, \dots, (w_1 - e_{\min})/\delta_e \end{cases} \Rightarrow C = 281 \times 141 = 39621 \quad (9) $$

where a manufacturability constraint was taken into account for the second level width ($w_2 \le w_1$). Each of these geometries are simulated with the FMM solver over the same spectral bandwidth as previously defined and for a TE state of polarization, which allows the 2 terms of the corresponding Taylor expansion to be represented in a phase space [22], as shown in Fig. 2(d). The color code represents the value of the mean transmission coefficient $|t|$ over the simulated bandwidth, when both heights $h_i$ are equal to 1 µm. With such a nano-element library, dispersion engineering consists in finding, for each position $x$ along the 1D metasurface, the shape of the nano-ridge that provides the best approximation of the requested phase function $\Phi_{\mathcal{F}}(x, \lambda)$. In a general way, this phase function can be also approximated by its Taylor expansion to the first order, namely

$$ \Phi_{\mathcal{F}}(x, \lambda) \approx \Phi_{\mathcal{F}}(x, \lambda_d) + (\lambda - \lambda_d) \left. \frac{\partial \Phi_{\mathcal{F}}}{\partial \lambda} \right|_{\lambda_d} \quad (10) $$

where the two coefficients $\Phi_{\mathcal{F}}(x, \lambda_d)$ and $\left. \frac{\partial \Phi_{\mathcal{F}}}{\partial \lambda} \right|_{\lambda_d}$ are function of the geometric characteristics of the problem, such as the pupil size $D_x$ of the 1D metasurface, the position ($f_x, f_z$) of the focal plane and its orientation $\alpha$ with respect to the metasurface (see Fig. 1).

To go further, it is necessary to define numerical values for these geometrical parameters, and as a first example, let us consider the case corresponding to a linear dispersion $a$ of 4 µm nm$^{-1}$, a pupil size $D_x$ of 175 µm, an axial distance $f_z$ equal to 14 times this pupil size ($f_z = 2.45$ mm), an offset distance $f_x$ equal to zero, and an orientation of the focal plane parallel to the metasurface ($\alpha = 0^\circ$). This orientation is the one which spontaneously comes to mind, so as to maximize the flux density reaching the focal plane. If we plot the two coefficients of the polynomial approximation (10) in the phase space of the nano-element library (see Fig. 2(e), green dots), we can see that the variation range of the specified phase derivative is roughly $\pm 0.6$ rad nm$^{-1}$, to be compared to the one covered by the library, i.e. $[-0.090, -0.065]$ rad nm$^{-1}$. This means that this library can not fulfill the phase specifications corresponding to this particular geometry ($f_x = 0, \alpha = 0^\circ$). This problem arises from the relatively low intrinsic dispersion of such nano-elements, which is a well known limitation to the design of dispersion engineered metasurfaces [20,22,26].

Therefore, to identify a configuration where, for each phase value $\Phi_{\mathcal{F}}(x, \lambda)$ defined by the functional $\mathcal{F}$, there is in the library a nano-element whose phase property $\phi(\lambda)$ correspond to this specified value with a very good approximation, we need to compute the two first coefficients of Taylor expansion of $\Phi_{\mathcal{F}}(x, \lambda)$, and verify that the corresponding point in the phase space is located in the dense part of the nano-element library. This dense part corresponds to the range of derivative values $(\partial\phi/\partial\lambda)_{\lambda_d}$ for which the condition $|\phi - \phi(\lambda_d)|<\epsilon$ is fulfilled, when $\phi$ is continuously varying between $-\pi$ and $+\pi$, and where $\epsilon$ is a very small phase quantity specified a
priori, e.g. $\pi/100$. As an example, in Fig. 2(d), this dense area corresponds to the derivative values between $-0.072 \text{ rad nm}^{-1}$ and $-0.084 \text{ rad nm}^{-1}$, i.e. a range $\Delta_{ne}$ of $0.012 \text{ rad nm}^{-1}$s. The use of a two-level nano-element library is a key point of this study. This allows the phase space to be filled in a sufficiently dense way so that the condition $|\varphi - \phi(\lambda_d)|<\pi/100$ is everywhere satisfied. Each level of the two nano-elements system is characterized by identical heights $h_1 = h_2 = 1 \text{ µm}$ and two independent variable widths $[w_1, w_2]$, whereas the one level counterpart has a height $h = h_1 + h_2 = 2 \text{ µm}$ and a single variable width $w$. Figure 2(c) shows the phase space filled with this one-level library, and it is clear that for such configuration, the filling of the phase space is not dense enough to satisfy the design condition for all points on the metasurface.

Due to the presence of a free term $C_{\theta,\lambda}$ in the expression of the functional $\mathcal{F}$ [see (1)], the constraint to satisfy becomes

$$ \Delta = \max \left( \left. \frac{\partial \Phi_{\mathcal{F}}(x, \lambda)}{\partial \lambda} \right|_{\lambda_d} \right) - \min \left( \left. \frac{\partial \Phi_{\mathcal{F}}(x, \lambda)}{\partial \lambda} \right|_{\lambda_d} \right) \le \Delta_{ne} \quad (11) $$

To illustrate this approach, consider the case where the offset distance $f_x$ is equal to its optimal value, i.e. $3.5 D_x$; we can calculate the variation of the parameter $\Delta$ versus the axial distance $f_z$ and the focal plane orientation angle $\alpha$. The result is shown in Fig. 3(a) (false color mapping). When $\alpha$ is greater than 75°, it becomes possible to design a metasurface over a range of axial distance $f_z$ that is a function of the orientation angle. To maximize the numerical aperture of the system, we need to minimize the ratio $f_z/D_x$, and a quite simple solution is provided by $\alpha = 90°$ and $f_z = 14 D_x$, i.e. a numerical aperture of 0.035. Figure 3(a) illustrates the importance of choosing an appropriate tilted angle in order to satisfy the inequality (11). For the given value of parameter $f_x$, the smaller the tilted angle (i.e when the detector tends to be parallel to the metasurface), the more dispersive the optical function $\Phi_{\mathcal{F}}$, and the more difficult it is to implement with our nano-element library. For $f_z$ up to $120 D_x$, it is impossible to satisfy (11) when $\alpha$ is smaller than 76°. Thus, the next results will always consider a detector with a tilted angle $\alpha = 90°$ which satisfy Eq. (11), a convenient choice that has already been chosen by [28]. Note that this approach confirms the impossibility to realize a spectral imaging function with a detector parallel to the metasurface ($\alpha = 0°$).

<figure><img src="image_6.png" /><figcaption>Fig. 3. Mapping of optical architecture that can be addressed by the nano-elements library. a) False color representation of the variation of parameter $\Delta$ as a function of the detector orientation $\alpha$ and the ratio between the axial distance $f_z$ and the pupil size $D_x$ ($D_x = 175 \text{ µm}$) - b) simplified layout of the optical system; the metasurface (MS) focuses and disperses a continuum of wavelengths on a detector (D) tilted by an angle $\alpha$ in the $(x, z)$ plane; the coordinates of the center of the detector are $(f_x, f_z)$ in the $(x, z)$ plane.</figcaption></figure>
4. Radiometry and point spread function

In Section 3, we have shown how a comparison between the spectral properties of two-level nano-elements stored in a library and the spectral dependence of the phase profile corresponding to the functional $\mathcal{F}$ allows the geometric characteristics of the hyperspectral imager to be optimized. However, knowing the resulting phase function $\Phi_{\mathcal{F}}(x, \lambda)$ does not provide any information about either the spectral dependence of the point spread function (PSF) of this imager along the focal plane nor on the irradiance level recorded at the top of these PSFs. The latter are quantities obtained from diffraction integrals and not only stigmatism condition from Fermat principle. To answer this question, we consider the system depicted on Fig. 4(a), where $\vec{u}_1$ (respectively, $\vec{u}_2$) is a unit vector perpendicular to the metasurface (respectively, to the focal plane), and $\rho$ is the distance between a fictitious emitter at position $x_e$ along the metasurface and a receiving point $(x, z)$ along the focal plane. For this 1D configuration and a TE polarization, the spatial distribution of the field amplitude along the focal plane corresponding to the ideal phase profile $\Phi_{\mathcal{F}}(x, \lambda)$ is given by the following Sommerfeld integral [37]:

$$ \mathcal{E}_y(x, z, \lambda) = \frac{i\bar{z}}{\sqrt{\lambda}} \int \frac{\exp\{i[k\rho + \phi_{\mathcal{F}}(x_e, \lambda)]\}}{\rho^{3/2}} dx_e \quad \text{where} \quad \rho = \sqrt{(x - x_e)^2 + z^2} \quad (12) $$

<figure><img src="image_7.png" /><figcaption>Fig. 4. a) Schematic representation of the optical system; the metasurface has an entrance pupil $D_x$ in the $(x, z)$ plane, and focus wavelengths on a matrix of photodiodes tilted by an angle $\alpha$ equal to 90°; from a ray tracing point of view, each ray from the metasurface is separated by a distance $\rho$ to the center of the detector, of coordinates $(f_x, f_z)$ - b) Plot of the ideal phase profile at the central wavelength - c) Illustration of a 2D off-axis zone plate designed to focus the central wavelength at the center of the detector. The middle slice of this 2D component has been chosen to define a 1D amplitude profile for the simulation (purple dashed line) - d) Absolute 1D irradiance calculated in the detector plane (D) considering the metasurface from $\lambda_{min}$ to $\lambda_{max}$ - e) Absolute 1D irradiance calculated in the detector plane (D) considering the ideal phase mask of b) - f) Absolute 1D irradiance calculated in the detector plane (D) considering the 1D zone plate of c).</figcaption></figure>
The flux density corresponding to this ideal case is obtained from the scalar product of the Poynting vector $\vec{\Pi}$ with the unit vector perpendicular to the detector plane $\vec{u}_2$, i.e:

$$ \frac{d\Phi}{dl} = \vec{\Pi} \cdot \vec{u}_2 = \Pi_x \sin \alpha + \Pi_z \cos \alpha = \frac{1}{2} \mathcal{E}_y^*(\mathcal{H}_z \sin \alpha - \mathcal{H}_x \cos \alpha) \quad (13) $$

$\mathcal{H}_x$ and $\mathcal{H}_z$ fields are calculated by differentiating Eq. (12) with respect to $z$ and $x$, respectively. Indeed, for a TE polarization:

$$ \mathcal{H}_x = -i\omega\mu_0 \frac{\partial \mathcal{E}_y}{\partial z} \quad ; \quad \mathcal{H}_z = i\omega\mu_0 \frac{\partial \mathcal{E}_y}{\partial x} \quad (14) $$

Fig. 4(b) illustrates the value of $\phi_{\mathcal{F}}(x_e, \lambda_d)$ along the fictitious surface at the central wavelength. Such parameter is used to define and ideal phase mask, which serves as a reference to compare the PSF obtained with the metasurface. We also define a 1D amplitude mask from an off-axis zone plate, focusing the central wavelength at the center of the detector (Fig. 4(c)).

## 5. Optical properties

In accordance with the conclusions of Section 3, the following geometric configuration is considered:

$$ D_x = 175 \ \mu\text{m} ; f_x = 3.5 D_x ; f_z = 14 D_x ; \alpha = 90^{\circ} ; \theta = 0^{\circ} \\ \lambda = [600, 700] \ \text{nm} ; a = 4 \ \mu\text{m} \text{nm}^{-1} $$

and the library of two-level nano-ridges presented in Fig. 2(a) is used to approximate the phase function of the corresponding spectral imager with a unit cell period of 400 nm.

### 5.1. Spectral response

The far field response of the 1D metasurface is calculated using a FMM solver with a number of orders equal to 2417, determined through a convergence test. Figure 4(d) shows the spatial variation of the 1D-irradiance along the u-axis obtained through this numerical calculation for 21 wavelengths regularly spaced from $\lambda_{min}$ to $\lambda_{max}$. For the same set of wavelengths, the irradiance is computed from a 1D Sommerfeld integral by considering the ideal phase mask (Fig. 4(e)) and the zone plate (Fig. 4(f)). The former use the quantities defined in Eq. (12) while the latter does not contain any phase inside the propagator term ($\phi_{\mathcal{F}}(x_e, \lambda) = 0$ for the zone plate) but rather a spatial amplitude term which is set accordingly into the integral. In each case, the intensity peak shifts from left to right as the wavelength decreases, as a natural consequence of a chromatic diffraction effect [20].

As expected, the 1D zone plate shapes the incident wavefront into a converging one, but the peak irradiances of the corresponding PSFs are far lower than those obtain with the ideal phase mask or the metasurface (intensity peak around 0.75 mW $\mu\text{m}^{-1}$, see Fig. 4(f)). Such decrease in this peak intensity is inherent to a zone plate in which a significant portion of the incoming light is either absorbed or reflected by the component (amplitude mask instead of phase mask). As an example, Fig. 5(a) compares the PSF of the designed metasurface (red curve) at the central wavelength $\lambda_d = 650$ nm with that obtained considering the ideal phase mask (blue curve). This small decrease in the peak intensity is due to a slight difference between the phase profile generated by the crossing of the metasurface and the analytic one defined by (7). The phase profiles generated by the metasurface at three wavelengths (600 nm, 650 nm, and 700 nm) are shown in Fig. 5(b), with a shape of tilted hyperbole in accordance with the off-axis nature of the optical function. Figures 5(c), 5(d), and 5(e) shows the phase difference between these profiles and the ideal function $\Phi_{\mathcal{F}}(\lambda)$ at the same wavelength. These slight differences justify a posteriori the dispersion engineering approach presented in Section 3.

In addition, the spatial full-width at half maximum FWHM of the PSF at 650 nm (34 $\mu\text{m}$) can be used to estimate the spectral resolution $\delta\lambda$ of the spectro-imager taking into account a linear dispersion $a$ of 4 $\mu\text{m} \text{nm}^{-1}$, i.e. $\delta\lambda = 8.5$ nm.
<figure><img src="image_8.png" /><figcaption>Fig. 5. a) Comparison between PSFs calculated at $\lambda_d = 650$ nm using 1D Sommerfeld integral with ideal phase mask (blue curve) and Fourier Modal Method (FMM) (red curve) - b) Phase profiles generated by the metasurface at three different wavelengths (600 nm, 650 nm, and 700 nm) - c) Residual difference between the phase profile generated by the metasurface at 600 nm and the ideal profile corresponding to the functional $\mathcal{F}$ - d) Same as graph c for $\lambda = 650$ nm - e) Same as graph c for $\lambda = 700$ nm.</figcaption></figure>

## 5.2. Spatial response

To determine the spatial response of the metasurface, it is necessary to consider different illumination angles $\theta$ (see Fig. 1(a)) and to analyze the evolution of the PSF (top irradiance and spatial width) along the v axis. However, since the metasurface used in our modeling being one-dimensional, the electromagnetic problem cannot be studied along y axis, which prevents having information on the v-dependence of the PSF. To nevertheless estimate the response of this metasurface in 2D and confirm that this component fulfills all the constraints corresponding to a hyperspectral pusbroom imager, we follow an approach directly derived from that used in Section 4 to analyze its spectral response. It is based on the use of a 2D Sommerfeld integral:

$$ \mathcal{E}_y(x, y, z, \theta, \lambda) = \frac{i}{\lambda} \iint \frac{|t|(x_e, y_e, \theta, \lambda) \exp\{i[k\rho + \Phi'(x_e, y_e, \theta, \lambda)]\}}{\rho} \cos(\vec{u}_1, \vec{\rho}) dx_e dy_e \quad (15) $$

where

$$ \rho = \sqrt{(x - x_e)^2 + (y - y_e)^2 + z^2} \quad \text{and} \quad \Phi'(x_e, y_e, \theta, \lambda) = \phi(x_e, y_e, \theta, \lambda) - ky_e \sin \theta \quad (16) $$
The function $\Phi'$ gathers the local nano-element response $\phi(x_e, y_e, \theta, \lambda)$ obtained with the FMM solver under different illumination angle $\theta$, and a phase shift term $-ky_e \sin\theta$ modeling the incoming tilted wavefront when considering a plane wave. As this theoretical component is two dimensional, the studied nano-elements are chosen to be a two-level squared shape pillars possessing the same geometrical parameters as the 1D counterpart (Fig. 6(a)). In a manner similar to what has been done before, each nano-element is represented in the phase space based on its phase response and dispersion (color points Fig. 6(b)). For a hyperspectral imaging system with a 2D component depicted on Fig. 6(c), the optical function at normal incidence $\Phi'(x_e, y_e, 0, \lambda)$ is Taylor expanded up to the first order, and plotted in the phase space (green dots Fig. 6(b)). This allows to choose a set of nano-elements to define the metasurface of aperture $D_x = D_y = 175 \ \mu m$. Each individual nano-element of the resulting metasurface is studied under different illumination angle $\theta$ along the bandwidth to get the global phase $\phi(x_e, y_e, \theta, \lambda)$ as well as the local transmission coefficient $|t|(x_e, y_e, \theta, \lambda)$ (Fig. 6(d)). The 2D PSF generated by the metasurface can then be estimated in the $(u, v)$ plane using a numerical computation of integral (15). Figure 7(a) shows the normalized irradiance map obtained for 5 illumination angles (0°, 2°, 4°, 6° and 8°). Irradiances are normalized with respect to the maximum peak intensity value of PSFs from 0° to 8°. We note the presence of a "smile" effect, i.e. a geometrical distortion of the PSF position for the same wavelength and different field angles: this effect is frequent in hyperspectral images. To cancel the contribution of this geometrical distortion effect, the evolution of the PSF peak intensity with respect to the illumination angle is plotted in Fig. 7(b) along a curve passing by the peak intensity at each illumination angle (magenta curve in Fig. 7(a)). Figure 6(c) and 6(d) show 1D graphs corresponding respectively to the normalized PSF variations with respect to the peak intensity, along $u$ axis for $v = 0$ and along $v$ axis for $u = 0$. From the 2D integral, a $10 \ \mu m$ FWHM for the PSF along the $v$ axis would be expected from this hypothetical 2D metasurface. From Fig. 7(a), the PSF shifts by an amount $dv = 100 \ \mu m$ with respect to a variation $d\theta = 2°$, leading to a spatial

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** a) schematic of the 2D nano-element building block under study for the design of a 2D metasurface. The geometrical parameters are identical to those of the 1D counterpart (TE polarisation) - b) comparative plots of the phase requirements of the metasurface (green dots) and the phase space addressed by the nano-element library of two-level square pillars [same color code as in Fig. 2(c)] - c) illustration of a 2D component with aperture ($D_x, D_y$) achieving the hyperspectral imaging function defined in Section 2 - d) pictures of phase maps $\phi(x_e, y_e, \theta, \lambda)$ and transmission maps $|t|(x_e, y_e, \theta, \lambda)$ used to approximate the far field response of the 2D metasurface based on Eq. (15).</figcaption></figure>
variation of the PSF $d\theta/dv = 0.02$ degree µm$^{-1}$. Taking into account the FWHM of the PSF along the *v* axis (FWHM = 10 µm) which limits the spatial resolution, such ideal 2D metasurface would have an angular resolution $\delta\theta = \frac{d\theta}{dv} \times \text{FWHM} = 0.2^\circ$.

<figure><img src="image_10.png" /><figcaption>Fig. 7. a) Irradiance map in the (u, v) plane showing the presence of a smile effect in the hyperspectral images - b) Normalized irradiance with respect to peak intensity at normal incidence $\theta = 0^\circ$ for different illumination angle up to 8° along the curvilinear smile axis - c) Normalized irradiance at normal incidence along the u axis - d) Normalized irradiance at normal incidence along the v axis.</figcaption></figure>

## 6. Conclusion

This work demonstrates the capability of a metasurface to simultaneously separate and focus spectral features in the scope of a pushbroom hyperspectral imager. This has been achieved through first, a full wave simulation restricted to a 1D problem, and second a confirmation and an extension of these modeling results by the calculation of 1D and 2D Sommerfeld integrals. The dispersion engineering approach can easily be extended to any spectral and angular optical function, provided that the functional $\mathcal{F}$ is properly defined with respect to $\theta$ and $\lambda$ dependencies. As for the case of achromatic lenses [22], the wavefront shaping at several wavelengths is inherently limited by the dispersion amplitude $\Delta_{ne}$ of the two-level nano-elements, restricting both numerical aperture and linear dispersion. To tackle this issue, the height of each level can be increased leading to a broader wavelength dispersion $d\phi/d\lambda$. Indeed, the nano-element library is implicitly based on an underlying propagation phase delay mechanism [38,39], in which the phase dispersion is proportional to the structure height [36]. However, increasing these heights leads to more complexity in the manufacturing process. As an alternative, it is possible to play with local resonances of meta-atoms [25] which increases the nano-element dispersion while keeping a realistic design to be fabricated. The dispersion engineering being validated for 1D metasurface, similar methodology could also apply to design metasurfaces components based on a library of 2D nano-elements composed of a two-level structure.

Finally, this work proposes a design methodology based on the relationship between properties at subwavelength scale and effective optical performances of the whole metasurface component. Such approach might offer to optical designer community new design pathways, providing easier understanding of metasurface performance for their introduction in future imaging architectures.
**Acknowledgments.** The authors acknowledge ANRT (Association Nationale de la Recherche et de la Technologie) for partial support of this work and their Thales colleagues Quentin Levesque, Thierry Viard and Hervé Benard for fruitful discussion about hyperspectral imaging systems.

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

1. N. V. Tkachenko, *Optical spectroscopy: methods and instrumentations* (Elsevier, 2006).

2. S. Mahesh, D. Jayas, J. Paliwal, and N. White, “Hyperspectral imaging to classify and monitor quality of agricultural materials,” *J. Stored Prod. Res.* **61**, 17–26 (2015).

3. W. Wang and J. Paliwal, “Near-infrared spectroscopy and imaging in food quality and safety,” *Sens. & Instrumen. Food Qual.* **1**(4), 193–207 (2007).

4. A. Y. Zhu, F. Yi, J. C. Reed, H. Zhu, and E. Cubukcu, “Optoelectromechanical multimodal biosensor with graphene active region,” *Nano Lett.* **14**(10), 5641–5649 (2014).

5. G. A. Shaw and H. K. Burke, “Spectral imaging for remote sensing,” *Linc. Lab. J.* **14**, 3–28 (2003).

6. F. Lemarquis, L. Abel-Tiberini, and C. Koc, “400-1000 nm all-dielectric linear variable filters for ultra compact spectrometers,” in *International Conference on Space Optics-ICSO 2010*, vol. 10565 (International Society for Optics and Photonics, 2019), p. 105655U.

7. K. Hendrix, “Linear variable filters for nasa’s ovirs instrument: pushing the envelope of blocking,” *Appl. Opt.* **56**(4), C201–C205 (2017).

8. H.-T. Chen, A. J. Taylor, and N. Yu, “A review of metasurfaces: physics and applications,” *Rep. Prog. Phys.* **79**(7), 076401 (2016).

9. Y. Zhao and A. Alù, “Manipulating light polarization with ultrathin plasmonic metasurfaces,” *Phys. Rev. B* **84**(20), 205428 (2011).

10. Y. Yang, W. Wang, P. Moitra, I. I. Kravchenko, D. P. Briggs, and J. Valentine, “Dielectric meta-reflectarray for broadband linear polarization conversion and optical vortex generation,” *Nano Lett.* **14**(3), 1394–1399 (2014).

11. M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, and F. Capasso, “Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging,” *Science* **352**(6290), 1190–1194 (2016).

12. X. Chen, L. Huang, H. Mühlenbernd, G. Li, B. Bai, Q. Tan, G. Jin, C.-W. Qiu, S. Zhang, and T. Zentgraf, “Dual-polarity plasmonic metalens for visible light,” *Nat. Commun.* **3**(1), 1198 (2012).

13. A. Arbabi, Y. Horie, M. Bagheri, and A. Faraon, “Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission,” *Nat. Nanotechnol.* **10**(11), 937–943 (2015).

14. F. Aieta, P. Genevet, M. A. Kats, N. Yu, R. Blanchard, Z. Gaburro, and F. Capasso, “Aberration-free ultrathin flat lenses and axicons at telecom wavelengths based on plasmonic metasurfaces,” *Nano Lett.* **12**(9), 4932–4936 (2012).

15. A. S. Glassner, *An introduction to ray tracing* (Elsevier, 1989).

16. S. D. Mellin and G. P. Nordin, “Limits of scalar diffraction theory and an iterative angular spectrum algorithm for finite aperture diffractive optical element design,” *Opt. Express* **8**(13), 705–722 (2001).

17. X. Li, S. Xiao, B. Cai, Q. He, T. J. Cui, and L. Zhou, “Flat metasurfaces to focus electromagnetic waves in reflection geometry,” *Opt. Lett.* **37**(23), 4940–4942 (2012).

18. D. A. Pommet, M. Moharam, and E. B. Grann, “Limits of scalar diffraction theory for diffractive phase elements,” *J. Opt. Soc. Am. A* **11**(6), 1827–1834 (1994).

19. S. J. Byrnes, A. Lenef, F. Aieta, and F. Capasso, “Designing large, high-efficiency, high-numerical-aperture, transmissive meta-lenses for visible light,” *Opt. Express* **24**(5), 5110–5124 (2016).

20. E. Arbabi, A. Arbabi, S. M. Kamali, Y. Horie, and A. Faraon, “Controlling the sign of chromatic dispersion in diffractive optics with dielectric metasurfaces,” *Optica* **4**(6), 625–632 (2017).

21. M. Khorasaninejad, F. Aieta, P. Kanhaiya, M. A. Kats, P. Genevet, D. Rousso, and F. Capasso, “Achromatic metasurface lens at telecommunication wavelengths,” *Nano Lett.* **15**(8), 5358–5362 (2015).

22. S. Shrestha, A. C. Overvig, M. Lu, A. Stein, and N. Yu, “Broadband achromatic dielectric metalenses,” *Light: Sci. Appl.* **7**(1), 85 (2018).

23. F. Aieta, M. A. Kats, P. Genevet, and F. Capasso, “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” *Science* **347**(6228), 1342–1345 (2015).

24. R. Sawant, P. Bhumkar, A. Y. Zhu, P. Ni, F. Capasso, and P. Genevet, “Mitigating chromatic dispersion with hybrid optical metasurfaces,” *Adv. Mater.* **31**(3), 1805555 (2019).

25. S. Wang, P. C. Wu, V.-C. Su, Y.-C. Lai, M.-K. Chen, H. Y. Kuo, B. H. Chen, Y. H. Chen, T.-T. Huang, and J.-H. Wang, “A broadband achromatic metalens in the visible,” *Nat. Nanotechnol.* **13**(3), 227–232 (2018).

26. W. T. Chen, A. Y. Zhu, J. Sisler, Z. Bharwani, and F. Capasso, “A broadband achromatic polarization-insensitive metalens consisting of anisotropic nanostructures,” *Nat. Commun.* **10**, 355 (2019).

27. M. Khorasaninejad, A. Y. Zhu, C. Roques-Carmes, W. T. Chen, J. Oh, I. Mishra, R. C. Devlin, and F. Capasso, “Polarization-insensitive metalenses at visible wavelengths,” *Nano Lett.* **16**(11), 7229–7234 (2016).
28. A. Y. Zhu, W. T. Chen, J. Sisler, K. M. Yousef, E. Lee, Y.-W. Huang, C.-W. Qiu, and F. Capasso, "Compact aberration-corrected spectrometers in the visible using dispersion-tailored metasurfaces," *Adv. Opt. Mater.* **7**(14), 1801144 (2019).

29. M. Faraji-Dana, E. Arbabi, A. Arbabi, S. M. Kamali, H. Kwon, and A. Faraon, "Compact folded metasurface spectrometer," *Nat. Commun.* **9**(1), 4196–4198 (2018).

30. R. Pestourie, C. Pérez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. G. Johnson, "Inverse design of large-area metasurfaces," *Opt. Express* **26**(26), 33732–33747 (2018).

31. J. S. Jensen and O. Sigmund, "Topology optimization for nano-photonics," *Laser Photonics Rev.* **5**(2), 308–321 (2011).

32. K. Svanberg, "A class of globally convergent optimization methods based on conservative convex separable approximations," *SIAM J. Optim.* **12**(2), 555–573 (2002).

33. M. Mansouree, H. Kwon, E. Arbabi, A. McClung, A. Faraon, and A. Arbabi, "Multifunctional 2.5 D metastructures enabled by adjoint optimization," *Optica* **7**(1), 77–84 (2020).

34. L. Li, "New formulation of the Fourier modal method for crossed surface-relief gratings," *J. Opt. Soc. Am. A* **14**(10), 2758–2767 (1997).

35. M. Moharam, E. B. Grann, D. A. Pommet, and T. Gaylord, "Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings," *J. Opt. Soc. Am. A* **12**(5), 1068–1076 (1995).

36. W. T. Chen, A. Y. Zhu, and F. Capasso, "Flat optics with dispersion-engineered metasurfaces," *Nat. Rev. Mater.* **5**(8), 604–620 (2020).

37. V. Nascov and P. C. Logofătu, "Fast computation algorithm for the Rayleigh-Sommerfeld diffraction formula using a type of scaled convolution," *Appl. Opt.* **48**(22), 4310–4319 (2009).

38. P. Lalanne and D. Lemercier-Lalanne, "On the effective medium theory of subwavelength periodic structures," *J. Mod. Opt.* **43**(10), 2063–2085 (1996).

39. P. Lalanne and P. Chavel, "Metalenses at visible wavelengths: past, present, perspectives," *Laser Photonics Rev.* **11**(3), 1600295 (2017).