Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Computational inverse design for cascaded systems of metasurface optics: comment

Sanaz Zarei and Amin Khavasi

Open Access Open Access

Abstract

In a recently published article by Backer [Opt. Express 27(21), 30308 (2019). [CrossRef]  ], a computational inverse design method is developed for designing optical systems composed of multiple metasurfaces. The forward propagation model used in this method was a discretized version of the angular spectrum propagator described by Goodman [Introduction to Fourier Optics, 1996]. However, slight modifications are necessary to increase the accuracy of this inverse design method. This comment examines the accuracy of the results obtained by the above-mentioned method by a full-wave electromagnetic solver and explains the reason of their difference. Thereafter, slight modifications to the method proposed by Backer are suggested, and the accuracy of final formulation is verified by a full-wave electromagnetic solver.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Backer combined the FDTD parameter sweep with Fourier-optics simulations to optimize a complex optical system composed of multiple metasurfaces with respect to individual meta-atom geometries [1]. This provided a relatively simple means of designing metasurface optics, without requiring rigorous FDTD simulations over the entire design volume. For his Fourier-optics simulations, he used a discretized version of the angular spectrum wave propagator described in [2]. Briefly the sampled electric field at the output of the m'th metasurface plane of the optical system is decomposed into a superposition of plane-wave components using the discretized Fourier transform, and then each plane wave component is multiplied by an appropriate phase factor that accounts for the phase accumulated as the wave propagates to the next surface within the optical system, and an inverse Fourier transform operation is then performed to determine the resulting electric field (Eq. (1)):

$${E^{out,{\lambda _j}}} = \left( {\prod\limits_{m = 1}^M {{F^\dagger }} {P^{m,{\lambda_j}}}F{\Phi ^{{w^m},{\lambda_j}}}} \right){E^{in,{\lambda _j}}}$$
which ${E^{in,{\lambda _j}}}$ is an (S×1)-dimensional vector containing discrete samples of the continuous input electric field distribution (at the wavelength ${\lambda _j}$), and each metasurface consists of S unique meta-atoms and there are total of M metasurfaces, ${E^{out,{\lambda _j}}}$ is a (S×1)-dimensional vector containing the resulting electric field at the output plane, and F and ${F^\dagger }$ are the discrete Fourier transform matrix and its inverse [1]. ${P^{m,{\lambda _j}}}$ is a diagonal matrix that effects the plane-wave propagation of the individual Fourier components of the incident electric field and ${\Phi ^{{w^m},{\lambda _j}}} = {e^{j{\phi ^{{w^m},{\lambda _j}}}}}$ is a diagonal matrix that contains complex exponentials associated with phase delays $({{\phi^{{w^m},{\lambda_j}}}} )$ induced by the choice of design variables ${w^m}$ at the input wavelength [1].

For optimization, a cost function in terms of the squared errors between a desired set of output intensity distributions, and realized set of intensity distributions is defined by using this model and is minimized iteratively by adjusting the design parameters $\{{{w^1},{w^2},\ldots ,{w^m}} \}$ using the steepest-descent algorithm [1]. In order to perform the optimization, the adjoint method in which the gradient of the cost function with respect to all design variables can be computed using only two full-field simulations, is utilized.

In this comment, the accuracy of the above-mentioned design technique is examined by a full-wave electromagnetic solver and the cause of the difference between their obtained results will be discussed and a few modifications will be made in the formulation to mitigate the error. Finally, the accuracy of the design method will be further increased by taking into account the transmission intensity of metasurfaces.

2. Verification with a full-wave electromagnetic solver and mitigation of error

To examine the accuracy of this method, a metasystem design example is chosen and its simulation results are compared with full-wave electromagnetic simulations. The chosen metasystem is composed of multi-layered cascaded 1D-metasurfaces. The 1D-metasurfaces (metalines) are 1D high-contrast transmitarray (1D HCTA) metasurfaces proposed in [3], which are composed of 1D rectangular etched slot arrays defined in a silicon-on-insulator substrate with silicon device layer thickness of 250nm and SiO2 insulator layer thickness of 2µm (Fig. 1(a)). The look-up table is generated by performing a parameter sweep over etched slots of different lengths, a constant slot width of 140nm, and a meta-atom pitch (a) of 500nm using commercial software package Lumerical FDTD, and measuring the output phase and intensity of the meta-atom for TE-polarized guided waves for the design wavelength of 1.55µm (Fig. 1(b)). To calculate the transmission phase and intensity of each meta-atom, FDTD ports were utilized with a fixed distance of 10µm along the x-direction.

 figure: Fig. 1.

Fig. 1. (a) Schematic of 1D rectangular etched slot arrays defined on the silicon-on-insulator (SOI) substrate, with a constant slot width (b) electric field transmission phase and amplitude of TE-polarized guided waves at the operation wavelength of 1.55 µm for a meta-atom versus slot length, while the slot width is fixed at 140 nm.

Download Full Size | PDF

A metasystem consisting of multilayered metalines in the SOI platform can be assumed as a two-dimensional problem. For a better verification, Lumerical 2D FDTD is used for parameter sweep and entire metasystem modelings. For all 2D simulations, the effective refractive index of the silicon slab waveguide ${n_{eff}}$ is used.

As a design example, a two-layered metalens composed of 200 meta-atoms on each layer (100µm aperture) will be designed using the computational inverse design proposed in [1]. The distance between the two layers is 200µm and it is desired that normally incident light comes to a focus of 200µm after the second metaline. The target intensity distribution (${I^{des,{\lambda _j}}}$) is selected as the ideal, diffraction-limited focal spot associated with a singlet metalens (one-layer metalens) with hyperboloid phase profile and focal length of 200µm and aperture size of 100µm. In the optimization, all the slot lengths within meta-atoms are initially set to 2µm.

Figure 2(a) shows the resulting optimized meta-atom lengths for both metalines, Fig. 2(b) shows the simulated intensity along the optical axis of metalens (along the x-axis) and Fig. 2(c) shows the simulated intensity at the focal length of metalens (along the y-axis). For verification, the designed doublet metalens is simulated using Lumerical 2D FDTD, which their obtained results are also shown in Fig. 2(b) and Fig. 2(c).

 figure: Fig. 2.

Fig. 2. (a) The optimized slot-lengths of first and second metalines for a metalens doublet design having 100 µm aperture. The simulated intensity that is obtained by different simulation methods (b) along the optical axis of metalens (along the x-axis) and (c) at the designed focal length of the metalens (along the y-axis). The small difference in the location of intensity peak between the Lumerical 2D FDTD and the other two methods in (b) is mainly due to the physical thickness of metalines in this software.

Download Full Size | PDF

As is clear, there is a large difference between the peak intensities of both methods. To get a better view of the problem, the designed metalens doublet is also simulated with the method utilized in [4,5]. This method uses full-wave simulations to obtain the local field for each meta-atom, and then assemble the field along the plane right after the metasurface using local periodic approximation and thereafter uses near-to-far-field transformation to calculate far-field distribution [5]. In this method, the distance between adjacent layers should be large enough so that only the far-field from one layer of the metasurface reaches the next layer. In both methods [1] and [5], if the fields are undersampled, then significant computational artifacts will result. The method proposed in [5], in contrary to the method proposed in [1] which is based on scalar wave approximation and therefore does not capture polarization effects, is not based on scalar diffraction theory and therefore it takes into account the polarization effects. However, the scalar-wave propagator proposed in [1] can be augmented to account for polarization by separately propagating the different polarization components, and incoherently summing the field intensities from the different polarizations to calculate the total intensity distribution at a given plane. The results obtained by method [5] is also depicted in Figs. 2(b) and 2(c). It is seen that these later results are better matched to 2D FDTD simulation results than the obtained results by the inverse design method presented in [1].

Our investigations indicate as the size of designed metasystem grows larger, i.e. the number of meta-atoms increases or the number of layers increases, the difference between the results obtained by the method proposed in [1] and the results achieved by Lumerical 2D FDTD increases. However, like the above-mentioned design, Lumerical results are more consistent with the analysis method proposed in [5].

To investigate the origin of the problem, we worked on the angular spectrum wave propagator described in [2]. The authors believe that in the angular spectrum wave propagator (Eq. (2) for 2D propagation), the integrals should be calculated from minus infinity to plus infinity and we are not eligible to restrict these integrals both in spatial domain and Fourier domain. In the case that we want to limit the spatial domain to only the metalines lengths, then the integral in Fourier domain should be calculated from minus infinity to plus infinity.

$$\tilde{E}(0,{f_y}) = \int\limits_{ - \infty }^\infty {E(0,y){e^{ - j2\pi {f_y}y}}dy} $$
$$E(x,y) = \int\limits_{ - \infty }^\infty {\tilde{E}(0,{f_y}){e^{j\frac{{2\pi }}{\lambda }\sqrt {1 - {{(\lambda {f_y})}^2}} x}}{e^{j2\pi {f_y}y}}d{f_y}} $$

To illustrate this issue, we assume that the transmission through a metaline is specified as $T(y) = A(y){e^{j\phi (y)}}$. Then limiting the Fourier integral of Eq. (2a) to only the metaline length is equivalent to multiplying the transmission function of the metaline to a rectangular function ($rect(y/L)$, where L is the length of the metaline):

$${T^{total}}(y) = T(y)rect(y/L)$$

Then the transmitted electric field at the plane right after the metaline can be expressed as:

$${E^T}(0,y) = T(y)rect(y/L){E^{in}}(0,y)$$
where ${E^{in}}(0,y)$ is the incident electric field to the metaline. By applying the convolution theorem, we have:
$${\tilde{E}^T}(0,{f_y}) = \tilde{T}({f_y}) \ast L\sin c(L \cdot {f_y})\ast {\tilde{E}^{in}}(0,{f_y})$$

As is known the sinc function expands from minus infinity to plus infinity and therefore for the most accurate results, the angular spectrum wave propagator integral in Eq. (2b) should be calculated from minus infinity to plus infinity. If we want to truncate the integral, then a suitable finite interval for ${f_y}$ should be selected.

In our case which the structure is quasi-periodic, sufficient number of Fourier samples are necessary for accurate results. Especially for larger x and y values in Eq. (2b), higher number of ${f_y}$ samples are required to achieve high accuracy. This means that for modeling larger metasystems (with larger number of neurons and/or larger distances between the metaline layers), higher number of Fourier samples are required.

However, in the inverse design procedure presented in [1], fast Fourier transform algorithm is used to simulate propagation of electric and adjoint fields due to its significant computational savings ($O(MS\log (S))$ operations for each forward propagation as opposed to $O(M{S^2})$ operations for an arbitrary series of M matrix multiplications [1]) and when using fast Fourier transform, the number of Fourier components is equal to the number of spatial samplings of the electric fields (which is equal to the number of meta-atoms (S)) and this is not sufficient to get accurate results. Therefore, while keeping the spatial resolution as is stated in [1] (equal to the number of meta-atoms (S)) in the simulations, the number of Fourier components should be increased further for accurate results. This means that after discretizing the angular spectrum wave propagator (Eq. (2b)), the Fourier sampling should be done as $f_{y}=\left(-\frac{1}{2 \times a}: \Delta f_{y}: \frac{1}{2 \times a}-\Delta f_{y}\right)$, where $\Delta {f_y} = \frac{1}{{{N_F} \times a}}$, a is the pitch of the metasurfaces (metalines) and NF is the number of Fourier components. Then, the accuracy can be enhanced by using a suitable value for NF (${N_F} > > S$).

We redesign the doublet metalens with this correction and verify the obtained results with Lumerical 2D FDTD and the method presented in [5]. We increase the number of Fourier components from 200 in previous simulation to 1200 in this simulation. The results are shown in Fig. 3. As is seen, now the results obtained with this method are closely matched with the method of [5] and is in good agreement with the results obtained by Lumerical 2D FDTD. However, this accuracy is achieved at the cost of lower calculation speed. In order to take the advantage of fast Fourier transform for fast calculations in metasystem design while having more accurate results, it is possible to enclose the metasurfaces with opaque screens (which their transmission is zero) and increase the number of meta-atoms far beyond the ones in metasurfaces (meta-atoms with zero transmission). This will dramatically increase the number of meta-atoms and therefore the number of Fourier components. This is equivalent to zero padding in digital signal processing. Zero-padding do the same function as increasing the number of Fourier samples, as zero-padding in one domain corresponds to a higher interpolation-density in the other domain [6]. In other words, you can interpolate the discrete Fourier transform (DFT) by zero padding, where zero padding enables you to obtain more accurate amplitude estimates of resolvable signal components [7].

 figure: Fig. 3.

Fig. 3. (a) The optimized slot-lengths of first and second metalines for a metalens doublet design having 100 µm aperture. The simulated intensity that is obtained by different simulation methods (b) along the optical axis of metalens (along the x-axis) and (c) at the designed focal length of the metalens (along the y-axis). The small difference in the location of intensity peak between the Lumerical 2D FDTD and the other methods in (b) is mainly due to the physical thickness of metalines in this software.

Download Full Size | PDF

The remaining difference between the results by two approximate analytical methods in [1] and [5] with Lumerical 2D FDTD (in Fig. 3) can be dominantly due to local periodic approximation used in these two methods, and secondarily due to the assumption of uniform transmission amplitude for all metasurface geometries and multiple Fresnel-like reflections of the intermediate layers of the optical system.

It should be mentioned that when a metasystem composed of multi-layered metasurfaces is forward propagated using the model proposed in [1], the results are almost accurate. The problem arises when a metasystem composed of multiple metasurfaces is designed using the inverse design method proposed in [1]. This is because there is a slight difference between the forward propagation results using the method proposed in [1] and the improved method presented in this comment. This slight difference doesn't make so much difference in the forward propagation of structures, but make a significant difference in the error back-propagation, which results in inaccurate inverse designed structures.

To give a better sense of the sufficient number of Fourier samples for the design procedure, a number of one-layer, two-layer and three-layer metalenses with 200 meta-atoms on each layer are inverse-designed using the method proposed in [1]. The propagation distance between layers is set to be equal to the focal length of metalenses and varies between 100µm and 900µm. The designed structures are forward propagated using different number of Fourier samples and the minimum required number of Fourier samples that results in the maximum error less than a threshold of 10−4 through the whole intensity profile is calculated. The minimum required number of Fourier samples at x1 = focus of the lens, y1 = middle point of metalines (50µm) and y2 = 5.25µm (which is very near to the edge of metalines) are computed (Fig. 4). As is seen in Fig. 4, the minimum required Fourier samples for one-layer, two-layer and three-layer structures are close to each other and mainly depends on the propagation distance. It can be inferred from these figures that for the metalens designs shown in Figs. 2 and 3, minimum of 1200 Fourier samples are required to have maximum error of 10−4 for the electric field intensity profile through the whole structure.

 figure: Fig. 4.

Fig. 4. (a) A two-layer metalens and the lines along which the minimum required number of Fourier samples are calculated. Minimum required number of Fourier samples for one-layer, two-layer and three-layer metalenses as a function of propagation distance are shown at (b) x1 = focus of the lens, (c) y1 = middle point of metalines (50 µm), and (d) y2 = 5.25 µm (nearly the edge of metalines).

Download Full Size | PDF

It is worth mentioning that some metasystems based on a more subwavelength meta-atom (a meta-atom pitch (a) of 350 nm), are also designed using the method proposed in [1]. We choose the constant slot width of 50 nm in these metasystems in order to the electric field transmission phase and amplitude of meta-atom versus slot length at 1.55 µm become very similar to the ones shown in Fig. 1(b). Similarly, for these metasystems, the difference between the results of method in [1] and other methods is significant. However, by decreasing the pitch from 500 nm to 350 nm, a general trend in the simulation results has not been observed.

As for the chosen architecture in this comment, the transmission amplitude is very close to 1 (Fig. 1(b)), the assumption of uniform transmission amplitude is a good approximation. However, in many metasurfaces, this may not be accurate. Although many metasurfaces in the literature are designed based on only the transmission phase of meta-atoms [8], better performing designs can be achieved by taking into account the transmission amplitude. We believe that the method proposed by Backer is a very efficient means for designing a system of multi-layered metasurfaces and to complete his formulation, in the next section, a few formulations are added to his one to take into account the transmission amplitude.

3. Increasing the accuracy of the model by assuming non-uniform transmission intensity for different meta-atom geometries

If we assume that any change in meta-atom geometry, not only affect the output phase, but also the output intensity, the forward propagation model will be modified as:

$${E^{out,{\lambda _j}}} = \left( {\prod\limits_{m = 1}^M {{F^\dagger }} {P^{m,{\lambda_j}}}F{T^{{w^m},{\lambda_j}}}{\Phi ^{{w^m},{\lambda_j}}}} \right){E^{in,{\lambda _j}}}$$
in which ${T^{{w^m},{\lambda _j}}}$ is a diagonal matrix containing the transmission amplitudes of meta-atoms on layer m. ${T^{{w^m},{\lambda _j}}}$ can be extracted from a look-up table like Fig. 1(b). For calculating $\frac{{dC}}{{d{w^m}}}$, the following expression can be obtained using the chain rule:
$$\frac{{dC}}{{d{w^m}}} = \frac{{dC}}{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \phi } }^{{w^m},{\lambda _j}}}}} \otimes \frac{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \phi } }^{{w^m},{\lambda _j}}}}}{{d{w^m}}} + \frac{{dC}}{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} }^{{w^m},{\lambda _j}}}}} \otimes \frac{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} }^{{w^m},{\lambda _j}}}}}{{d{w^m}}}$$
where ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} ^{{w^m},{\lambda _j}}}$ is an (S×1)-dimensional vector containing the diagonal elements of ${T^{{w^m},{\lambda _j}}}$ matrix. The symbol ${\otimes} $ denotes the element-wise multiplication. The first expression on the right hand side of Eq. (7) can be calculated from the procedure presented in [1] while including ${T^{{w^m},{\lambda _j}}}$ in the formulations. For calculating the second expression, $\frac{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} }^{{w^m},{\lambda _j}}}}}{{d{w^m}}}$ can be easily obtained from Fig. 1(b), while for computing $\frac{{dC}}{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} }^{{w^m},{\lambda _j}}}}}$ at a single spatial location s’ on the m'th metasurface, we use the chain rule again:
$$\frac{{dC}}{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} }_{s^{\prime}}}^{{w^m},{\lambda _j}}}} = \sum\limits_{s = 1}^S {\left( {\frac{{dC}}{{dE_s^{out,{\lambda_j}}}}} \right)} \left( {\frac{{dE_s^{out,{\lambda_j}}}}{{d\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T}_{s^{\prime}}^{{w^m},{\lambda_j}}}}} \right) = 4R\left\{ {{{({{E^{out,{\lambda_j}}} \otimes ({I^{{\lambda_j}}} - {I^{des,{\lambda_j}}})} )}^\dagger }\left( {\frac{{dE_{}^{out,{\lambda_j}}}}{{d\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T}_{s^{\prime}}^{{w^m},{\lambda_j}}}}} \right)} \right\}$$
Where Writinger Calculus for computing derivatives with respect to a complex variable is employed. $R\{{} \}$ denotes the real part of the expression enclosed in brackets, ${I^{{\lambda _j}}}$ and ${I^{des,{\lambda _j}}}$ are realized and desired intensity distributions, and $^\dagger $ indicates the adjoint of a matrix. The derivative $\frac{{dE_{}^{out,{\lambda _j}}}}{{d\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} _{s^{\prime}}^{{w^m},{\lambda _j}}}}$ can be calculated using Eq. (6) as:
$$\frac{{dE_{}^{out,{\lambda _j}}}}{{d\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} _{s^{\prime}}^{{w^m},{\lambda _j}}}} = \left( {\prod\limits_{m^{\prime} = m + 1}^M {{F^\dagger }{P^{m^{\prime},{\lambda_j}}}F{T^{{w^{m^{\prime}}},{\lambda_j}}}{\Phi ^{{w^{m^{\prime}}},{\lambda_j}}}} } \right)({{F^\dagger }{P^{m,{\lambda_j}}}F{\Phi ^{{w^m},{\lambda_j}}}} )\Delta _{^2}^{s^{\prime}}{E^{m - 1,{\lambda _j}}}$$
in which ${E^{m - 1,{\lambda _j}}}$ is the intermediate incident electric field upon the m'th metasurface and $\Delta _{^2}^{s^{\prime}}$ is a diagonal matrix containing only one non-zero element:
$$\Delta _{2,\xi ,\xi }^{s^{\prime}} = \left\{ {\begin{array}{{l}} {1\;\quad \quad \quad if\;\xi = s^{\prime}}\\ {0\quad \quad \quad otherwise} \end{array}} \right.$$

By substituting Eq. (9) into Eq. (8) and taking the adjoint of the resulted equation, we have:

$$\frac{{dC}}{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} }^{{w^m},{\lambda _j}}}}} = 4R\left\{ {{{({E^{m - 1,{\lambda_j}}}^\dagger )}^T} \otimes ({{\Phi ^{{w^m},{\lambda_j}}}^\dagger {F^\dagger }{P^{m,{\lambda_j}}}^\dagger F} )\left( {\prod\limits_{m^{\prime} = 0}^{M - m - 1} {{\Phi ^{{w^{M - m^{\prime},{\lambda_j}}}}}^\dagger {T^{{w^{M - m^{\prime},{\lambda_j}}}}}^\dagger {F^\dagger }{P^{M - m^{\prime},{\lambda_j}}}^\dagger F} } \right){a^{{\lambda_j}}}} \right\}$$
where ${a^{{\lambda _j}}}$ is the adjoint field vector and is defined as
$${a^{{\lambda _j}}} = {E^{out,{\lambda _j}}} \otimes ({I^{{\lambda _j}}} - {I^{des,{\lambda _j}}})$$

By substituting Eq. (11) into Eq. (7), the gradient $\frac{{dC}}{{d{w^m}}}$ can be efficiently computed.

To show how this design framework can be more efficient in the design of metasurface systems than the one that only takes into account the transmission phase of meta-atoms, we design two metalens doublets based on the transmission phase and amplitude characteristics presented in Fig. 5(b). These transmission characteristics are achieved in a similar configuration to Fig. 1(a) (Fig. 5(a)), by performing a parameter sweep over etched slots of different widths, a constant slot length of 1.25µm, and a meta-atom pitch (a) of 500nm.

 figure: Fig. 5.

Fig. 5. (a) Schematic of 1D rectangular etched slot arrays defined on SOI substrate, with a constant slot length (b) electric field transmission phase and amplitude of TE-polarized guided waves at the operation wavelength of 1.55 µm for a meta-atom versus slot width, while the slot length is fixed at 1.25 µm.

Download Full Size | PDF

One metalens is designed by using the framework which only takes into account the phase of meta-atoms (Figs. 6(a)-(c)), and the other metalens is designed by using the presented design framework in this comment, which takes into account both transmission phase and amplitude of meta-atoms (Figs. 6(d)-(f)). The designed metalens doublets are composed of 200 meta-atoms on each layer (100µm aperture). The distance between the two layers is set to be 200µm and it is desired that normally incident light comes to a focus of 200µm after the second metaline (like previous designs). The target intensity distribution (${I^{des,{\lambda _j}}}$) is selected as the ideal, diffraction-limited focal spot associated with a singlet metalens (one-layer metalens) with hyperboloid phase profile and focal length of 200µm and aperture size of 100µm. In the optimizations, all the slot widths within meta-atoms are initially set to be 270nm and due to fabrication constraints, the minimum and maximum slot widths are specified as 50nm and 450nm, respectively.

 figure: Fig. 6.

Fig. 6. Designed metalens doublet having 100 µm aperture by only taking into account the transmission phase (P): (a) The optimized slot-widths of the first and second metalines. The simulated intensity that is obtained by different simulation methods (b) along the optical axis of metalens (along the x-axis) and (c) at the designed focal length of the metalens (along the y-axis). Designed metalens doublet having 100 µm aperture by taking into account both the transmission phase and amplitude (P&A): (d) The optimized slot-widths of the first and second metalines. The simulated intensity that is obtained by different simulation methods (e) along the optical axis of metalens (along the x-axis) and (f) at the designed focal length of the metalens (along the y-axis). The small difference in the location of intensity peak between the Lumerical 2D FDTD and the other methods in (b) and (e) is mainly due to the physical thickness of metalines in this software.

Download Full Size | PDF

As is seen in Fig. 6(b) and Fig. 6(c), for the metalens that is designed only based on transmission phase, the focus intensity drops significantly when simulated using both the transmission phase and amplitude information (compared to when simulated using only the transmission phase). However, as is shown in Figs. 6(e) and 6(f), the metalens that is designed based on both transmission phase and amplitude, shows higher focus intensity (and in general higher performance). For a comparison, we also simulate the metalens of Figs. 6(d)-(f) using only transmission phase information.

4. Conclusion

In conclusion, the shortcoming of the method presented in [1] is when using fast Fourier transform for computing the plane wave expansion of electric fields at the plane right after the metasurfaces, as the number of Fourier components (plane waves) in this case is equal to the number of meta-atoms and this is not sufficient for accurate results.

Funding

Iran National Science Foundation (4002443).

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. A. Backer, “Computational inverse design for cascaded systems of metasurface optics,” Opt. Express 27(21), 30308 (2019). [CrossRef]  

2. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996).

3. Z. Wang, T. Li, A. Soman, D. Mao, T. Kananen, and T. Gu, “On-chip wavefront shaping with dielectric metasurface,” Nat. Commun. 10(1), 3547 (2019). [CrossRef]  

4. Z. Wu, M. Zhou, E. Khoram, B. Liu, and Z. Yu, “Neuromorphic metasurface,” Photonics Res. 8(1), 46 (2020). [CrossRef]  

5. 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 (2018). [CrossRef]  

6. J. O. Smith, “Zero Padding Applications,” in Mathematics of the Discrete Fourier Transform (DFT) with Audio Applications, 2nd Edition, http://ccrma.stanford.edu/∼jos/mdft/Fourier_Theorems_DFT.html, online book, 2007 edition, accessed <14/3/2022 > .

7. https://www.mathworks.com/help/signal/ug/amplitude-estimation-and-zero-padding.html.

8. N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater. 13(2), 139–150 (2014). [CrossRef]  

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.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (6)
Fig. 1.
Fig. 1. (a) Schematic of 1D rectangular etched slot arrays defined on the silicon-on-insulator (SOI) substrate, with a constant slot width (b) electric field transmission phase and amplitude of TE-polarized guided waves at the operation wavelength of 1.55 µm for a meta-atom versus slot length, while the slot width is fixed at 140 nm.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. (a) The optimized slot-lengths of first and second metalines for a metalens doublet design having 100 µm aperture. The simulated intensity that is obtained by different simulation methods (b) along the optical axis of metalens (along the x-axis) and (c) at the designed focal length of the metalens (along the y-axis). The small difference in the location of intensity peak between the Lumerical 2D FDTD and the other two methods in (b) is mainly due to the physical thickness of metalines in this software.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. (a) The optimized slot-lengths of first and second metalines for a metalens doublet design having 100 µm aperture. The simulated intensity that is obtained by different simulation methods (b) along the optical axis of metalens (along the x-axis) and (c) at the designed focal length of the metalens (along the y-axis). The small difference in the location of intensity peak between the Lumerical 2D FDTD and the other methods in (b) is mainly due to the physical thickness of metalines in this software.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. (a) A two-layer metalens and the lines along which the minimum required number of Fourier samples are calculated. Minimum required number of Fourier samples for one-layer, two-layer and three-layer metalenses as a function of propagation distance are shown at (b) x1 = focus of the lens, (c) y1 = middle point of metalines (50 µm), and (d) y2 = 5.25 µm (nearly the edge of metalines).

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. (a) Schematic of 1D rectangular etched slot arrays defined on SOI substrate, with a constant slot length (b) electric field transmission phase and amplitude of TE-polarized guided waves at the operation wavelength of 1.55 µm for a meta-atom versus slot width, while the slot length is fixed at 1.25 µm.

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. Designed metalens doublet having 100 µm aperture by only taking into account the transmission phase (P): (a) The optimized slot-widths of the first and second metalines. The simulated intensity that is obtained by different simulation methods (b) along the optical axis of metalens (along the x-axis) and (c) at the designed focal length of the metalens (along the y-axis). Designed metalens doublet having 100 µm aperture by taking into account both the transmission phase and amplitude (P&A): (d) The optimized slot-widths of the first and second metalines. The simulated intensity that is obtained by different simulation methods (e) along the optical axis of metalens (along the x-axis) and (f) at the designed focal length of the metalens (along the y-axis). The small difference in the location of intensity peak between the Lumerical 2D FDTD and the other methods in (b) and (e) is mainly due to the physical thickness of metalines in this software.

View in Article | Download Full Size | PDF

Equations (13)

Equations on this page are rendered with MathJax. Learn more.

E o u t , λ j = ( ∏ m = 1 M F † P m , λ j F Φ w m , λ j ) E i n , λ j
E ~ ( 0 , f y ) = ∫ − ∞ ∞ E ( 0 , y ) e − j 2 π f y y d y
E ( x , y ) = ∫ − ∞ ∞ E ~ ( 0 , f y ) e j 2 π λ 1 − ( λ f y ) 2 x e j 2 π f y y d f y
T t o t a l ( y ) = T ( y ) r e c t ( y / L )
E T ( 0 , y ) = T ( y ) r e c t ( y / L ) E i n ( 0 , y )
E ~ T ( 0 , f y ) = T ~ ( f y ) ∗ L sin ⁡ c ( L ⋅ f y ) ∗ E ~ i n ( 0 , f y )
E o u t , λ j = ( ∏ m = 1 M F † P m , λ j F T w m , λ j Φ w m , λ j ) E i n , λ j
d C d w m = d C d ϕ ⇀ w m , λ j ⊗ d ϕ ⇀ w m , λ j d w m + d C d T ⇀ w m , λ j ⊗ d T ⇀ w m , λ j d w m
d C d T ⇀ s ′ w m , λ j = ∑ s = 1 S ( d C d E s o u t , λ j ) ( d E s o u t , λ j d T ⇀ s ′ w m , λ j ) = 4 R { ( E o u t , λ j ⊗ ( I λ j − I d e s , λ j ) ) † ( d E o u t , λ j d T ⇀ s ′ w m , λ j ) }
d E o u t , λ j d T ⇀ s ′ w m , λ j = ( ∏ m ′ = m + 1 M F † P m ′ , λ j F T w m ′ , λ j Φ w m ′ , λ j ) ( F † P m , λ j F Φ w m , λ j ) Δ 2 s ′ E m − 1 , λ j
Δ 2 , ξ , ξ s ′ = { 1 i f ξ = s ′ 0 o t h e r w i s e
d C d T ⇀ w m , λ j = 4 R { ( E m − 1 , λ j † ) T ⊗ ( Φ w m , λ j † F † P m , λ j † F ) ( ∏ m ′ = 0 M − m − 1 Φ w M − m ′ , λ j † T w M − m ′ , λ j † F † P M − m ′ , λ j † F ) a λ j }
a λ j = E o u t , λ j ⊗ ( I λ j − I d e s , λ j )
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use