Simple strategy for the simulation of axially symmetric large-area metasurfaces

Augusto Martins; Augusto Martins; Achiles F. da Mota; Chris Stanford; Taylor Contreras; Justo Martin-Albo; Alexander Kish; Carlos O. Escobar; Adam Para; Roxanne Guenette

doi:10.1364/JOSAB.510725

1. INTRODUCTION

Metalenses consist of subwavelength nanostructures that locally modify the phase profile of an incoming beam to focus it [1–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 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–12]. However, the larger the metalens area, the more computationally expensive it is to perform accurate simulations. 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 times the operating wavelength [16]. However, the most common approach to simulate even larger area metalenses is the local approximation method [6,14], 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. Inspiration can be drawn from the simulation of diffractive optical elements (DOEs), which is already an established technique. Common approaches involve performing zone decomposition to account for edge effects [17,18] and local grating approximation [19]. Although accurate for DOE simulation, these models could not fully simulate a metalens given its local control of the light phase and polarization.

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

Download Full Size | PDF

Here, we propose a method to accurately simulate large-area metalenses in a timely manner at single wavelength operation. Our strategy is based on sampling the metalens phase gradient profile and modeling 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 polarization. A similar approach has been recently demonstrated to simulate quantum emitters’ response near periodic patterned hyperbolic metamaterials [20,21]. 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 used our model to simulate metalenses made of 1200 nm tall fused silica circular nanoposts placed 500 nm apart with different diameters used to control the phase profile operating at a wavelength of 632 nm. We compared our method against rigorous FDTD simulations and managed to reduce the simulation time and memory requirements by at least one order of magnitude while keeping the results accurate. Section 3.A of the main text contains an analysis of a small metalens with an aperture diameter of 100 µm, and relevant figures of merit are obtained and compared against FDTD simulations. In Section 3.B we use the model to calculate the local diffraction efficiency of a metalens as a function of position and angle of incidence with 100 µm wide Gaussian beam excitation. The metalens has a focal length of 7 mm and numerical aperture (NA) of 0.68, and the model is again compared against FDTD simulations. In Section 3.C, we show the simulation results of a metalens with 2.26 mm of diameter using the wave-optics and ray tracing models, which would be impractical to simulate on FDTD. Finally, in Section 3 of Supplement 1, we even demonstrated the simulation of a doublet metalens using the ray tracing model.

2. LOCAL LINEAR GRATING DIFFRACTION ORDER DECOMPOSITION

We propose a model based on a transfer-function approach that considers coupling among different posts. Let us assume a metalens whose design phase profile is given by $\psi (\vec r)$, where $\vec r$ is the in-plane radial vector from its center. The phase gradient is sampled on $N$ positions, ${\vec r_a},a \in [1,N]$, as illustrated in Fig. 1. The region around the point ${\vec r_a}$ can be approximated by as a blazed binary grating with period ${P_a}$ given by the grating vector equation as

(1)$${P_a} = \frac{{2\pi}}{{G({{\vec r}_a})}};\quad G(\vec r) = \left| {\frac{{\partial \psi (\vec r)}}{{\partial r}}} \right|,$$

where $G({\vec r_a})$ is the radial phase gradient value calculated at ${\vec r_a}$ and is identical to the corresponding blazed binary grating wave vector ${G_a}$. 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).

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 [22,23]. The transfer function is defined as the diffracted field amplitudes of the transmitted/reflected waves as a function of the incoming angle of incidence, in spherical coordinates, written in terms of the in-plane wave-vector components’ amplitudes. We solve Maxwell’s equations for a supercell of the blazed binary grating on a region with a cross-section area given by $d \times {P_a}$, where $d$ 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 ${P_a}$ 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 following the corresponding metasurface nanostructure design. 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 accelerated using accurate and faster methods [14–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. On the wave-vector model, we use a finite number of approximations along the azimuthal direction. We address the impact of this approximation for increasing distances from the metalens vertex in Supplement 1. We summarized the method in the block diagram shown in Fig. 2. The diagram summarizes the linear grating approximation described in this section, and the wave-vector and ray tracing methods, which are described in Sections 2.A and 2.B, respectively.

Fig. 2. Block diagram summarizing the proposed method.

Download Full Size | PDF

Fig. 3. (a) shows how the different sectors are set up based on the gradient sampling. (b) shows the angular splitting for each sector, forming a given region. $\vec r$ is defined as the position vector centered on the metalens, and ${\vec r_{\textit{ab}}}$ is the vector pointing to the center of the region $ab$. ${\hat \rho _{\textit{ab}}}$ and ${\hat \phi _{\textit{ab}}}$ are the blazed direction and the corresponding orthogonal unit vectors of region $ab$, respectively.

Download Full Size | PDF

A. Vectorial Wave Implementation of the New 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 [24]. As explained in the previous section, we sample the metalens phase profile gradient with linear patches as shown in Fig. 3(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. 3(b) and 3(c). That is, we approximate the metalens phase profile by a piece-wise linear phase profile at the points ${\vec r_{a,b}}$ creating $I$ sectors radially, and we split each sector into ${J_a}$ regions where $a$ indicates a given radial sector, as represented in Fig. 3. The total transmitted or reflected fields by the metalens are calculated as the sum of the contribution of each sector

(2)$$\vec E(\vec r) = \sum\limits_{a = 1}^I \sum\limits_{b = 1}^{{J_a}} {\vec E_{\textit{ab}}}(\vec r){W_{\textit{ab}}}(\vec r),$$

where ${\vec E_{a,b}}(\vec r)$ is the electric field distribution calculated on region $b$ of sector $a$, $\vec r$ is defined as a radial vector centered on the metalens vertex [see Fig. 3(a)], and ${W_{\textit{ab}}}(\vec r)$ is a window function, as shown in Fig. 3(c), that limits the region area and is given by Eq. (3):

(3)$${W_{\textit{ab}}}(\vec r) = \left\{{\begin{array}{ll}{1,}&\;\;{{\rm if}\,\,{r_{a - 1}} \lt |\vec r| \le {r_a} \wedge \left| {\phi - \frac{{2\pi b}}{{{J_a}}}} \right| \le \frac{{2\pi}}{{{J_a}}}}\\{0,}&\;\;{{\rm otherwise}}\end{array}} \right..$$

If the incoming field distribution has a Fourier transform given by ${{\boldsymbol {\vec E}}_0}({\vec k_\parallel})$ (bold symbols represent the Fourier transformed quantities), the transmitted field angular spectrum on axial region $b$ of sector a (see Fig. 3) is given by

(4)$${{\boldsymbol {\vec E}}_{\textit{ab}}}({\vec k_\parallel}) = \int d{\kappa ^2}{\mathop T\limits^ \leftrightarrow _{\textit{ab}}}({\vec k_\parallel},{\vec \kappa _\parallel}) \cdot {{\boldsymbol{ \vec E}}_0}({\vec \kappa _\parallel}),$$

where ${\vec k_\parallel}$ and ${\vec \kappa _\parallel}$ are, respectively, the incoming field and dummy in-plane reciprocal vectors. Equation (4) 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 [20,21,24]. The tensor ${\mathop T\limits^ \leftrightarrow _{\textit{ab}}}({\vec k_\parallel},{\vec \kappa _\parallel})$ is the transfer function and is given by

(5)$$\begin{split}{\mathop T\limits^ \leftrightarrow _{\textit{ab}}}({\vec k_\parallel},{\vec \kappa _\parallel}) &= \sum\limits_g \sum\limits_{\alpha ,\beta} \delta ({\vec k_\parallel} - {\vec \kappa _\parallel} - {\vec G_{\textit{abg}}}){{\cal T}_{abg\alpha \beta}}({\vec \kappa _\parallel}){\hat p_\alpha}({\vec \kappa _\parallel} \\&\quad+ {\vec G_{\textit{abg}}}){\hat p_\beta}({\vec \kappa _\parallel}),\end{split}$$

where ${{\cal T}_{abg\alpha \beta}}({\vec \kappa _\parallel})$ is the transmission coefficient function of the $a$th blazed grating at the $b$th sector for the diffraction order $g$ and incident polarization ${{\hat p}_\beta}(\vec \kappa)$ and diffracted polarization ${{\hat p}_\alpha}({\vec \kappa _\parallel} + {\vec G_{\textit{abg}}})$. Note that the “$\parallel$” subindex indicates a component of a vector in the metalens plane. The dummy indices $\alpha$ and $\beta$ can assume either $s$ or $p$ for the polarization states, which are defined as

(6)$${\hat p_s}({\vec k_\parallel}) = \frac{{{{\vec k}_\parallel} + \hat z{k_ \bot}}}{{|{{\vec k}_\parallel} + \hat z{k_ \bot}|}} \times \hat z,$$

(7)$${\hat p_p}({\vec k_\parallel}) = \frac{{{{\vec k}_\parallel} + \hat z{k_ \bot}}}{{|{{\vec k}_\parallel} + \hat z{k_ \bot}|}} \times {\hat p_s}({\vec k_\parallel}).$$

$\hat z$ is normal to the metalens plane, and ${k_ \bot}$ is the normal component of the wave vector and can be found according to the dispersion equation for homogeneous media

(8)$${k_ \bot} = \sqrt {{n^2}k_0^2 - |{{\vec k}_\parallel}{|^2}} ,$$

where $n$ is the surrounding medium refractive index and ${k_0}$ the free space wavenumber. The rotation symmetry of the metalens phase profile allows us to perform the following approximation:

(9)$${\mathop T\limits^ \leftrightarrow _{\textit{ab}}}({\vec k_\parallel},{\vec \kappa _\parallel}) = {R_{{\phi _b} - {\phi _1}}}\left[{{{\mathop T\limits^ \leftrightarrow}_{a{\phi _1}}}({{\vec k}_\parallel},{{\vec \kappa}_\parallel})} \right].$$

We define ${R_\phi}$ as the rotation operator on the metalens plane by an angle $\phi$, and ${\phi _b}$ the central polar angle of region $b$, assuming that the blazed grating is aligned with region 1. The wave vector of the diffraction order $G$ is given by the grating equation

(10)$${\vec G_{\textit{abg}}} = {m_g}{G_{\textit{ar}}}{\hat \rho _{\textit{ab}}} + {n_g}{G_\phi}{\hat \phi _{\textit{ab}}},{m_g},{n_g} \in \mathbb{Z}.$$

${\hat \rho _{\textit{ab}}}$ and ${\hat \phi _{\textit{ab}}}$ are the radial and azimuthal vectors at the region $b$ in sector $a$ [see Fig. 2(c)] obtained at the point ${\vec r_{\textit{ab}}}$, where the approximation is made. Note that the dummy index $g$ maps all diffraction orders. ${G_{\textit{abg}}}$ and ${G_\phi} = \frac{{2\pi}}{d}$ ($d$ is the metasurface unit cell) are the blazed grating wave vectors along the blazing and orthogonal directions, respectively. ${{\cal T}_{abg\alpha \beta}}({\vec \kappa _\parallel})$ is the transmission coefficient matrix obtained rigorously from Maxwel’s equations using the RCWA method with plane wave excitation from the substrate side with incoming in-plane wave vector ${\vec \kappa _\parallel}$. The Dirac delta function $\delta ({\vec k_\parallel} - {\vec \kappa _\parallel} - {\vec G_{\textit{abg}}})$ models the additional linear phase profile in real space from the blazed grating diffraction (grating equation). Thus, from Eqs. (4) and (5), and after integrating over the Dirac delta function, the output field spectrum is given by

(11)$$\begin{split}{{{{\boldsymbol {\vec E}}}_{{\boldsymbol a\boldsymbol b}}}({{\vec k}_\parallel}) }&={ \sum\limits_g \sum\limits_{\alpha ,\beta} {{\cal T}_{abg\alpha \beta}}({{\vec k}_\parallel} - {{\vec G}_{\textit{abg}}}){{\hat p}_\alpha}({{\vec k}_\parallel}){{\hat p}_\beta}({{\vec k}_\parallel} - {{\vec G}_{\textit{abg}}})} \\&\quad\cdot{ {{{\boldsymbol {\vec E}}}_0}({{\vec k}_\parallel} - {{\vec G}_{\textit{abg}}})},\end{split}$$

where each term in the summation of Eq. (11), ${{\boldsymbol {\vec e}}_{\textit{abg}}}({\vec k_\parallel})$, is the spectrum of the output field produced by the $g$th diffraction order and is given by

(12)$${{\boldsymbol {\vec e}}_{\textit{abg}}}({\vec k_\parallel}) = \sum\limits_{\alpha ,\beta} {{\cal T}_{abg\alpha \beta}}({\vec k_\parallel}){\hat p_\alpha}({\vec k_\parallel} + {\vec G_{\textit{abg}}}){\hat p_\beta}({\vec k_\parallel}) \cdot {{\boldsymbol {\vec E}}_0}({\vec k_\parallel}).$$

Finally, to obtain the field in real space, we simply perform an inverse Fourier transform in ${{\boldsymbol {\vec e}}_{\textit{abg}}}({\vec k_\parallel})$. Applying the inverse Fourier transform and the shifting property in Eq. (11) we have that

(13)$${\vec E_{\textit{ab}}}({\vec r_\parallel}) = \sum\limits_g {e^{- i{{\vec G}_{\textit{abg}}} \cdot {{\vec r}_\parallel}}}\int dk_\parallel ^2{{\boldsymbol {\vec e}}_{\textit{abg}}}({\vec k_\parallel}){e^{i{{\vec k}_\parallel} \cdot {{\vec r}_\parallel}}}.$$

That is, the total field on the ($a$, $b$) region is given by the linear superposition of the diffracted fields, as one would expect, modulated by the linear phase of the order. Note that from Eq. (13), the phase term ${e^{- j{{\vec G}_{\textit{abg}}} \cdot {{\vec r}_\parallel}}}$ controls the central position of the diffracted field spectrum and creates a local spatial linear phase profile, 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 metasurface. To eliminate this problem, we correct this phase term by substituting it to the original phase profile, $\psi ({\vec r_\parallel})$, as follows:

(14)$${\vec E_{\textit{ab}}}({\vec r_\parallel}) = \sum\limits_g {e^{- i(- {m_g}\phi ({{\vec r}_\parallel}) + {n_g}G\phi |{{\vec r}_\parallel}|)}}\int dk_\parallel ^2{{\boldsymbol {\vec e}}_{\textit{abg}}}({\vec k_\parallel}){e^{i{{\vec k}_\parallel} \cdot {{\vec r}_\parallel}}}.$$

Note that the phase gradient of Eqs. (13) and (14) should be the same at the phase gradient approximation point ${\vec r_{\textit{ab}}}$. This approximation simply corrects the phase value to adjust to the original curvature of the phase profile. Finally, any deviation on the phase profile due to wavefront errors by the metalens is captured by the higher-order modes. This approximation is valid if the phase gradient sampling is sufficiently small. In our model, we sampled the gradient profile in points and obtained a good match between FDTD simulations and our method.

B. Ray Tracing Implementation of the New Model

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. Here we describe how to calculate the scattered rays by the metalens upon illumination with an arbitrary ray. Our main goal is to analyze the efficiency of the diffraction efficiency of the metalens using geometrical optics. This approach makes the calculation of the metalens efficiency much faster than the diffraction method described in the previous section. Furthermore, we envisage the integration of this method into other ray tracing compatible packages to allow for fast metalens simulations, such as GEANT4, which is a platform for the simulation of the passage of particles through matter using Monte Carlo methods and enables metalens integration on particle physics simulations [25]. The metalens is now treated as a phase discontinuity on the ray path [26] but also accounting for the diffraction efficiency of the metalens. Given a ray with in-plane wave vector ${\vec \kappa _\parallel}$ incoming at a point ${\vec r_\parallel}$ on the metalens [see Fig. 3(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 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 ${T_G}({\vec k_\parallel},\vec P)$ and ${R_G}({\vec k_\parallel},\vec P)$ as being the probabilities of a ray with incoming in-plane wave vector ${\vec k_\parallel}$ and polarization state $\vec P$ being diffracted in transmission and reflection, respectively, and $G \in \mathbb{D}$ is a given diffraction order in the set of stored diffraction orders space ($\mathbb{D}$). We describe in Supplement 1 how to find the probabilities from the RCWA simulations. From energy conservation, we can define the probability of a ray being absorbed as

(15)$$A({\vec \kappa _\parallel},\vec P) = 1 - \sum\limits_{G \in \mathbb{D}} {R_G}({\vec \kappa _\parallel},\vec P) + {T_G}({\vec \kappa _\parallel},\vec P).$$

Note that $P \equiv a\hat s({\vec \kappa _\parallel}) + b\hat p({\vec \kappa _\parallel})$ is the Jones vector on the linear polarization basis given by Eqs. (6) and (7). With the probabilities of each outcome (reflection, transmission or absorption) we define the following cumulative probability distribution function $f$ as

(16)$$f[i] = \left\{{\begin{array}{*{20}{l}}{\sum\nolimits_{j = 1}^i {V[j]}}&{1 \le i \le M + N + 1}\\0&{i = 0}\end{array}}, \right.$$

where $V[i] \equiv [{T_{{G_1}}},{T_{{G_2}}}, \cdots {T_{{G_N}}},{R_{{G_1}}},{R_{{G_2}}}, \cdots {R_{{G_M}}},A]$ and we omitted the region indices for clarity. Therefore, we can use $f$ to define intervals where a given outcome might happen. The diffraction order can then be found by generating a random number and analyzing in which interval it has fallen. Given the uniformly distributed random variable $\chi \in [0,1]$, we can obtain $i$, and consequently the outcome of the event from the sequence $V$, by solving the following equation:

(17)$$g(\chi) = i,\,\,{\rm if}\,\,f[i - 1] \lt \chi \le f[i],$$

where $g$ is the resulting diffraction index and index of the sequence $V$, which is used to map a given event. With the scattering event determined, the resulting ray can be either 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 wave vector given by

(18)$$\vec k = {\vec \kappa _\parallel} + {\vec G_a} \pm \hat n{k_n},$$

where $\hat n$ is the normal vector to the metasurface, ${\vec G_a}$ is the reciprocal vector of the $g$th order, and ${k_n}$ 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 wave vector. As shown in Supplement 1, we can also obtain the diffracted polarization by the metalens. Each ray can be labeled 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 Supplement 1, where we show the output ray order decomposition for a normal incident ray distribution on the metalens. Such discrimination is also possible on the vectorial wave method, but it would require storing a field distribution of each order, which is impractical.

3. RESULTS

A. Focusing Profile Comparison

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 $y = 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.C. We use 1200 nm tall fused silica glass-based nanoposts in our design [12,27]. The metalens focal length is 200 µm with a diameter of 100 µm (NA = 0.25), and it encodes a hyperbolic phase profile, which produces a perfect spherical output beam and can reach diffraction-limited focusing [4]. The vectorial wave based model qualitatively reproduced the FDTD field distribution with good fidelity, even accounting for the appearance of higher-order focal spots as shown in Figs. 4(a) and 4(b), respectively, at normal incidence. The model also reproduced the field distribution at ${20^\circ}$ of incidence as shown by Figs. 4(e) and 4(f), respectively. We calculated the normalized root mean square (RMS) of the field distribution intensity difference between our vectorial wave model and FDTD integrated over the $xz$ plane. We normalized the field distributions to the integral of the FDTD field distributions on the $XZ$ plane. The normalized field difference intensity RMS is $1.9 \times {10^{- 3}}$ and $5.2 \times {10^{- 3}}$ at normal and 20° of incidence. The ray tracing model also simulates the focusing field distribution and the high-order focal spots as shown in Figs. 4(c) and 4(g). However, it accounts only for the power flow of the diffracted orders and idealizes the intensity distribution as they are modeled as geometrical rays, which do not manifest interference effects. The interference of the diffraction orders 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. 4(d) and 4(h) at normal and oblique (${20^ \circ}$) incidence, respectively. These curves were normalized with to their own power flux and peak intensity of the ideal Airy disk distribution. The field distributions at the focal plane obtained by our method present 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. 4(h). We also calculated the focusing efficiency on a circle with radius three times larger than the ideal Airy disk full width at half maximum [28]. At normal incidence, we obtained 73.1%, 71.2%, and 71.5% with the FDTD, vectorial wave, and ray tracing methods, respectively. Finally, the ray tracing model allows us to easily tag each ray according to the diffraction order it originated from as shown in Figs. S1–S3 of Supplement 1, where we simulate different metalens systems, including a doublet one.

Fig. 4. 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)–(c) [(e)–(g)] show the longitudinal field amplitude distributions focused by the metalens at normal [${20^ \circ}$] incidence calculated using the vectorial wave linear patch, FDTD, and ray tracing with linear patching. (d) and (h) show the transversal cuts of the field distributions on the focal plane at normal and ${20^ \circ}$ incidence, respectively.

Download Full Size | PDF

B. 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 fused silica glass-based nanopost design described in [27] 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 times the operating wavelength) waist. The simulation region for each scan point is $200\,\,\unicode{x00B5}\rm m\times 200\,\,\unicode{x00B5}\rm m$.

The simulation results at normal ${10^\circ}$ and ${20^\circ}$ of incidence are shown in Figs. 5(a)–5(e), Figs. 5(f)–5(j), and Figs. 5(k)–5(o), respectively. The calculated diffraction efficiencies of orders $m = - {2}$, ${-}{1}$, and 0; $m = {2}$; and $n = {0}$ correspond to each column of Fig. 5, where the metalens focusing comes from the $m = - {1}$ order. All methods agree for all simulated angles of incidence and diffraction orders. In particular, the efficiency of the $m = - 1$ order is the highest at normal incidence, and it is almost constant throughout the whole scan, remaining higher than 60%, which is a 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. Additional experimental data agrees well with the simulation; see [27].

Fig. 5. Diffraction efficiencies of different ($m$, $n$) 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 and ${10^ \circ}$ and ${20^ \circ}$ 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.

Download Full Size | PDF

C. Computational Resources

In this section, we compare the simulation time and memory of each method requirement against FDTD simulations to calculate the transmission through the metalens without free space propagation. 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 [27], 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 \times2R \times 2\,\,\unicode{x00B5}\rm{m} $, where ${ R}$ is the metalens radius $(R = f\tan ({\rm asin}(\rm NA))$. 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 [24]. The FDTD simulations are carried out on the FAS cluster with 128 cpus distributed over four 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 accelerate 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\,\,{\rm rays}/{{\unicode{x00B5}{\rm m}}^2}$ hitting the metalens. Figures 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.

Fig. 6. (a) and (b) show, respectively, the memory and time required to simulate fused silica glass-based metalenses with different radius ($R$) (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 four nodes using parallelization. The other models were run in serial and have room for improvement. All simulations are performed on a normally incident plane wave operating at 632 nm.

Download Full Size | PDF

4. LARGE-AREA METALENS FOCUSING PROFILE SIMULATIONS

The low computational time and minimal memory usage of our approach 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 ${\sim}3755{\lambda _0}$. The metalenses are excited by plane waves at normal and ${20^ \circ}$ of incidence. Figures 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 more than 100 h and require more than a petabyte of memory to simulate on a FDTD cluster with more than 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.

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^ \circ}$ of incidence, respectively. (c) and (d) [(g) and (h)] show the focusing profile of a hyperbolic metalens at normal and ${20^ \circ}$ of incidence, respectively. The metalenses’ focal lengths are 1 mm with a diameter of 2.26 mm (NA = 0.75), and the operating wavelength is 632 nm.

Download Full Size | PDF

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}\;{\unicode{x00B5}{\rm m}}$, 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 [29].

5. 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 using either a ray tracing approach or a vectorial wave model, depending on the desired application. We found very good agreement between our model and FDTD simulations. This approach reduces the time and memory requirements by orders of magnitude, allowing the simulation of metalenses with diameters larger than $3755{\lambda _0}$ on a personal computer.

Funding

U.S. Department of Energy (L2021.011, DE-AC02-07CH11359); 2021 Alfred P. Sloan Research Fellowship (ERC-2020-SyG-951281).

Acknowledgment

Most computations in this paper were run on the FASRC cluster supported by the FAS Division of Science Research Computing Group at Harvard University. This document was prepared by using in part the resources of the Fermi National Accelerator Laboratory (Fermilab) and the Noble Liquid Test Facility (NLTF), a U.S. Department of Energy, Office of Science, Office of High Energy Physics HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA). A. Martins and R. Guenette were partly supported by the Alfred P. Sloan Foundation.

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.

Supplemental document

See Supplement 1 for supporting content.

REFERENCES

1. P. Lalanne and P. Chavel, “Metalenses at visible wavelengths: past, present, perspectives,” Laser Photonics Rev. 11, 1600295 (2017). [CrossRef]

2. M. Khorasaninejad and F. Capasso, “Metalenses: Versatile multifunctional photonic components,” Science 358, eaam8100 (2017). [CrossRef]

3. H. Liang, A. Martins, B.-H. V. Borges, et al., “High performance metalenses: numerical aperture, aberrations, chromaticity, and trade-offs,” Optica 6, 1461–1470 (2019). [CrossRef]

4. M. Khorasaninejad, W. T. Chen, R. C. Devlin, et al., “Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging,” Science 352, 1190–1194 (2016). [CrossRef]

5. B. Groever, W. T. Chen, and F. Capasso, “Meta-lens doublet in the visible region,” Nano Lett. 17, 4902–4907 (2017). [CrossRef]

6. A. Martins, K. Li, J. Li, et al., “On metalenses with arbitrarily wide field of view,” ACS Photonics 7, 2073–2079 (2020). [CrossRef]

7. M. Khorasaninejad, Z. Shi, A. Y. Zhu, et al., “Achromatic metalens over 60 nm bandwidth in the visible and metalens with reverse chromatic dispersion,” Nano Lett. 17, 1819–1824 (2017). [CrossRef]

8. H. Pahlevaninezhad, M. Khorasaninejad, Y.-W. Huang, et al., “Nano-optic endoscope for high-resolution optical coherence tomography in vivo,” Nat. Photonics 12, 540–547 (2018). [CrossRef]

9. J.-S. Park, S. Zhang, A. She, et al., “All-glass, large metalens at visible wavelength using deep-ultraviolet projection lithography,” Nano Lett. 19, 8673–8682 (2019). [CrossRef]

10. J. Kim, J. Seong, W. Kim, et al., “Scalable manufacturing of high-index atomic layer-polymer hybrid metasurfaces for metaphotonics in the visible,” Nat. Mater. 22, 474–481 (2023). [CrossRef]

11. L. Zhang, S. Chang, X. Chen, et al., “High-efficiency, 80 mm aperture metalens telescope,” Nano Lett. 23, 51–57 (2023). [CrossRef]

12. J.-S. Park, S. W. D. Lim, A. Amirzhan, et al., “All-glass 100 mm diameter visible metalens for imaging the cosmos,” ACS Nano 18, 3187–3198 (2024). [CrossRef]

13. T. W. Hughes, M. Minkov, V. Liu, et al., “A perspective on the pathway toward full wave simulation of large area metalenses,” Appl. Phys. Lett. 119, 150502 (2021). [CrossRef]

14. R. Pestourie, C. Pérez-Arancibia, Z. Lin, et al., “Inverse design of large-area metasurfaces,” Opt. Express 26, 33732–33747 (2018). [CrossRef]

15. H.-C. Lin, Z. Wang, and C. W. Hsu, “Fast multi-source nanophotonic simulations using augmented partial factorization,” Nat. Comput. Sci. 2, 815–822 (2022). [CrossRef]

16. J. Skarda, R. Trivedi, L. Su, et al., “Low-overhead distribution strategy for simulation and optimization of large-area metasurfaces,” npj Comput. Mater. 8, 78 (2022). [CrossRef]

17. A. Nemes-Czopf, D. Bercsényi, and G. Erdei, “Simulation of relief-type diffractive lenses in zemax using parametric modelling and scalar diffraction,” Appl. Opt. 58, 8931–8942 (2019). [CrossRef]

18. A. Nemes-Czopf and G. Erdei, “Combined finite-element and scalar diffraction simulation of light scattering on zone edges of diffractive intraocular lenses,” Appl. Opt. 62, 6491–6498 (2023). [CrossRef]

19. C. Jang, O. Mercier, K. Bang, et al., “Design and fabrication of freeform holographic optical elements,” ACM Trans. Graph. 39, 184 (2020). [CrossRef]

20. A. F. da Mota, A. Martins, H. Ottevaere, et al., “Semianalytical model for design and analysis of grating-assisted radiation emission of quantum emitters in hyperbolic metamaterials,” ACS Photonics 5, 1951–1959 (2018). [CrossRef]

21. A. F. da Mota, A. Martins, V. Pepino, et al., “Semianalytical modeling of arbitrarily distributed quantum emitters embedded in nanopatterned hyperbolic metamaterials,” J. Opt. Soc. Am. B 36, 1273–1287 (2019). [CrossRef]

22. D. M. Whittaker and I. S. Culshaw, “Scattering-matrix treatment of patterned multilayer photonic structures,” Phys. Rev. B 60, 2610–2618 (1999). [CrossRef]

23. E. Popov and M. Nevière, “Maxwell equations in Fourier space: fast-converging formulation for diffraction by arbitrary shaped, periodic, anisotropic media,” J. Opt. Soc. Am. A 18, 2886–2894 (2001). [CrossRef]

24. L. Novotny and B. Hecht, Principles of Nano-optics (Cambridge University, 2006).

25. S. Agostinelli, J. Allison, K. A. Amako, et al., “Geant4-a simulation toolkit,” Nucl. Instrum. Methods Phys. Res. A 506, 250–303 (2003). [CrossRef]

26. N. Yu, P. Genevet, M. A. Kats, et al., “Light propagation with phase discontinuities: Generalized laws of reflection and refraction,” Science 334, 333–337 (2011). [CrossRef]

27. T. Contreras, A. Martins, C. Stanford, et al., “A method to characterize metalenses for light collection applications,” J. Instrum. 18, T09004 (2023). [CrossRef]

28. R. Menon and B. Sensale-Rodriguez, “Inconsistencies of metalens performance and comparison with conventional diffractive optics,” Nat. Photonics 17, 923–924 (2023). [CrossRef]

29. X. Ni, S. Ishii, A. V. Kildishev, et al., “Ultra-thin, planar, babinet-inverted plasmonic metalenses,” Light Sci. Appl. 2, e72 (2013). [CrossRef]

Simple strategy for the simulation of axially symmetric large-area metasurfaces

Abstract

1. INTRODUCTION

2. LOCAL LINEAR GRATING DIFFRACTION ORDER DECOMPOSITION

A. Vectorial Wave Implementation of the New Model

B. Ray Tracing Implementation of the New Model

3. RESULTS

A. Focusing Profile Comparison

B. Diffraction Efficiency Calculation

C. Computational Resources

4. LARGE-AREA METALENS FOCUSING PROFILE SIMULATIONS

5. CONCLUSION

Funding

Acknowledgment

Disclosures

Data availability

Supplemental document

REFERENCES

Supplementary Material (1)

Data availability

Cited By

Figures (7)

Equations (18)

Journal of the Optical Society of America B