K-series approximation of vectorial optical fields for designing diffractive optical elements with subwavelength feature sizes

TL;DR

This paper introduces a wave-vector series approximation method for vectorial optical fields in subwavelength diffractive optical elements, enabling efficient and accurate design by addressing computational challenges of vectorial simulations.

Contribution

It proposes a k-series approximation algorithm based on RCWA for efficient vectorial field simulation in subwavelength DOEs, outperforming scalar models in accuracy.

Findings

Predicted intensity profiles match FDTD simulations with fractional error.

Algorithm achieves 10-100 times better accuracy than scalar models.

Requires reasonable computational resources for complex designs.

Abstract

Diffractive optical elements (DOEs) are widely applied as compact solutions for desired light manipulations via wavefront shaping. Recent advanced chip applications further require their feature sizes to move down to the subwavelength, which inevitably brings forth vectorial effects of optical fields and makes the typical scalar-based theory invalid. However, simulating and optimizing their vectorial fields, which are associated with billions of adjustable parameters in the optical element, are difficult to do, because of the issues of numerical stability and the highly-demanding computational cost. To address this problem, this research proposes an applicable algorithm by means of a wave-vector (k) series approximation of vectorial optical fields. On the basis of the semi-analytical rigorous coupled wave analysis (RCWA), an adequate selection scheme on k-series enables computationally efficient yet still predictive calculations for DOEs. The performance estimations for exemplary designs by the finite difference time domain (FDTD) method show that the predicted intensity profiles by the proposed algorithm agree with the target by just a fractional error. Together with optimizing the geometrical degrees of freedom (e.g., DOE depth h) as compensation for errors from the truncation of k-series, the algorithm demonstrates its outperformance by one or two orders of magnitude in accuracy versus the scalar-based model, and demands only a reasonable computational resource.

Equations

$$\epsilon (\mathrm{x},\mathrm{y})$$

$${\epsilon }^{-1}(\mathrm{x},\mathrm{y})$$

$$\mathbf{E}$$

$$\mathbf{H}$$

$${\mathrm{k}}_{\mathrm{xu}}$$

$${\mathrm{k}}_{\mathrm{yv}}$$

$${\mathrm{k}}_{\mathrm{z},\mathrm{uv}}$$

$${\partial }_{\mathrm{z}}{\mathrm{e}}_{\mathrm{x},\mathrm{uv}}(\mathrm{z})$$

$${\partial }_{\mathrm{z}}{\mathrm{e}}_{\mathrm{y},\mathrm{uv}}(\mathrm{z})$$

$${\partial }_{\mathrm{z}}{\mathrm{h}}_{\mathrm{x},\mathrm{uv}}(\mathrm{z})$$

$${\partial }_{\mathrm{z}}{\mathrm{h}}_{\mathrm{y},\mathrm{uv}}(\mathrm{z})$$

$${\mathrm{R}}_{\mathrm{x}(\mathrm{y})}({\mathrm{x}}_{\mathrm{I}},{\mathrm{y}}_{\mathrm{I}})\subset \mathrm{F}{\left[{\mathrm{P}}_{\mathrm{x}(\mathrm{y})}({\mathrm{x}}_{\mathrm{D},\mathrm{s}},{\mathrm{y}}_{\mathrm{D},\mathrm{s}})\right]}^{-1}$$

$$\mathrm{I}({\mathrm{x}}_{\mathrm{I}},{\mathrm{y}}_{\mathrm{I}})={\left|\mathbf{R}({\mathrm{x}}_{\mathrm{I}},{\mathrm{y}}_{\mathrm{I}})\right|}^2.$$

$$\mathrm{m}=\int \int {\left|{\mathrm{I}}_0({\mathrm{x}}_{\mathrm{I}},{\mathrm{y}}_{\mathrm{I}})-\alpha \mathrm{I}({\mathrm{x}}_{\mathrm{I}},{\mathrm{y}}_{\mathrm{I}})\right|}^2\mathrm{d}{\mathrm{x}}_{\mathrm{I}}\mathrm{d}{\mathrm{y}}_{\mathrm{I}}.$$

$$\alpha =\frac{\int \int {\mathrm{I}}_0({\mathrm{x}}_{\mathrm{I}},{\mathrm{y}}_{\mathrm{I}})\mathrm{d}{\mathrm{x}}_{\mathrm{I}}\mathrm{d}{\mathrm{y}}_{\mathrm{I}}}{\int \int \mathrm{I}({\mathrm{x}}_{\mathrm{I}},{\mathrm{y}}_{\mathrm{I}})\mathrm{d}{\mathrm{x}}_{\mathrm{I}}\mathrm{d}{\mathrm{y}}_{\mathrm{I}}}.$$

$$\mathrm{\Theta }(\mathrm{\Delta m},\mathrm{T})={\mathrm{e}}^{-\frac{\mathrm{\Delta m}}{\mathrm{T}}},$$

$$\mathrm{T}=\frac{{\mathrm{T}}_0}{1+\mathrm{n}},1\le \mathrm{n}\le {\mathrm{N}}_{\mathrm{A}}.$$

Full text

††journal: oe††articletype: Research Article

1 Introduction

Diffractive optical elements (DOEs) are devices composed of transparent materials with complex topography to realize specified functions (e.g. optical computing and cell trapping [1, 2]) that are not feasible with standard refractive optics. With the rapid development of additive manufacturing techniques [3, 4, 5] in recent years, especially two-photon polymerization lithography (TPL) [6, 7, 8], it has become more feasible to fabricate DOEs with feature sizes down to ten nanometers. So far this technology is widely used for optical applications in the visible band, such as holograms for virtual reality [9], colorful 3D prints [10], and optical anticounterfeiting devices [11]. Moreover, TPL has been proven to bear potential at producing DOEs for advanced chip-scale applications like all-optical diffractive deep neural networks [12, 13, 14]. However, for the mentioned feature size, the assumptions of typical scalar-based design methods are violated [15, 16] and the vectorial effects of optical fields become more pronounced [17]. A general, fully rigorous design tool is therefore needed to enable an accurate synthesis of such elements. So far this is quite a challenging work for researchers.

A small strand of the literature has reported on designing DOEs with subwavelength feature sizes. For DOEs in the 1D domain, D.W. Prather et al. proposed an optimization-based algorithm using the boundary element method (BEM) [18]. Its computational costs, however, prevent the design of realistic DOEs in reasonable time frames. J. Jiang et al. developed the micro-genetic algorithm FDTD method [19], but it remains a challenging design tool due to its time-consuming searching routines. M.E. Testorf et al. demonstrated the Gerchberg-Saxton-algorithm FDTD method [20], which is a highly efficient optimization algorithm, but is related to the paraxial design algorithm. Feng Di et al. employed an iterative optimization algorithm [21, 22] based on the rigorous electromagnetic theory for a 1D beam splitter, but it is still subject to time-consuming and local-searching-only issues. In the 2D domain, H. Hao et al. presented a hybrid method that combines the rigorous coupled wave analysis (RCWA) [23, 24] and the genetic algorithm [25] to design 2D wide-angle beam splitters. This method, however, precludes designs from DOE’s applications other than splitters. S. J. Byrnes et al. demonstrated RCWA-based simulations to realize novel meta-lenses [26]. The technique takes advantage of the peripherally-radial periodicity and the circular symmetry of a lens to efficiently design the metasurfaces, but could treat other symmetries in an unfavorable way; e.g., C2 symmetry in binary DOEs or asymmetry in ternary DOEs. S. Wei et al. utilized optical properties of the unit cell (an elliptical silicon pillar in a $\mathrm{0.5\times 0.5\times 0.6\ \mu m^{3}}$ box) to design the dielectric metasurface of the polarimeter system [27], featuring high (above $\mathrm{60\%}$) detection efficiency for polarized light. This semi-analytical approach, however, could inherently suppress further improvement of the design due to strong phase gradients, interactions between neighboring pillars, and oblique angles [26]. N. Li et al. explored a large-area ($\mathrm{2500\mu m\times 2500\mu m}$) metasurface beam deflector [28], showing the pixel level ($\mathrm{120\mu m\times 120\mu m}$) controlling light. As discussed by the authors, the finite efficiency of the device needs to be improved by optimizing the design of the nanopillars, taking the interactions among nearby pillars with various diameters into consideration. Thus, it is an urgent requirement to develop an efficient algorithm for designs of DOEs that are finite in extent (mm- to cm-scale), have subwavelength features, and are aperiodic for different functions.

On the back of the scalar-based simulated annealing algorithm (s-SAA) [29], this research proposes a vector-based simulated annealing algorithm (v-SAA) for designing DOEs with subwavelength feature sizes. By considering the vectorial effects of optical fields dominatingly taking place inside the elements, this work mainly re-models the wave front function $\mathrm{U}$ in s-SAA with the rigorous electromagnetic theory. Specifically, on the basis of RCWA, a k-series approximation method is introduced to quickly calculate the pixel-field profile $\mathrm{\mathbf{P}}$ (in contrast to $\mathrm{U}$ in s-SAA) during the v-SAA process, where the order of k-series correlates with optional parameters for selecting different degrees of accuracy and computational costs. The v-SAA hence can achieve computationally efficient yet still predictive optical designs for DOEs.

The rest of the paper runs as follows. Four subsections of section 2 are arranged based on the following considerations: identify problems in the process, find solutions to problems, integrate the solutions into the process, and improve the new assembly (with new solution component) process. Specifically, section 2 describes the detailed design technique and points out the limits of the scalar-based theory on subwavelength-scale designs. To break the barriers of scalar-based inaccuracy and the vector-based high cost, a k-series approximation method (KSAM) is introduced to efficiently compute the vectorial optical fields of given DOE morphologies $\mathrm{\Upsilon}$. On the basis of KSAM, an iterative algorithm v-SAA is then realized for efficient designs of DOEs with subwavelength feature sizes. The complementary optimization process for errors from the truncation of k-series is finally mentioned for high-quality products. Section 3 presents the design results and analysis to demonstrate the effectiveness of the proposed method. Section 4 concludes this work.

2 Design principle

2.1 Limits of the scalar-based theory

To interpret the deviation from the prediction by the scalar-based theory in systems with subwavelength feature size, in this subsection we apply the rigorous electromagnetic theory, FDTD, for accurate computations of vectorial optical fields. The vectorial effects are then demonstrated by observing an exemplary DOE at different feature sizes $\mathrm{\Lambda_{ph}}$. The DOE is composed of $5\times 5$ square pixels with feature size $\mathrm{\Lambda_{ph}}$ and depth h, as indicated in Fig. 1(c). Here two kinds of media, the substrate medium with refractive index $n_{2}=1.5$ and air with $n_{1}=1.0$, occupy the pink and transparent pixels, respectively. Figure 1(a) illustrates the transmission fields on the exiting surface of DOE (at z=h+$\mathrm{\delta_{h}}$=$\mathrm{\lambda+\delta_{h}}$), providing a plane-wave incidence with wavelength $\lambda$ at $\mathrm{z=0}$ and structure parameter $\mathrm{\Lambda_{ph}/\lambda=6}$. An appropriate separation $\mathrm{\delta_{h}}$ from the surface $\mathrm{z=h}$ is adopted to diminish the influence of evanescent waves. In this classical regime of refractive optics- i.e., $\mathrm{\Lambda_{ph}/\lambda\gg 1}$ - numerical results show that the intensity profile of optical fields is only dependent on the morphology (the pink or transparent pixel) of an individual pixel and has almost constant amplitude value $\mathrm{|E|=1}$. In Fig. 1(a’) the corresponding phase profile $\mathrm{arg(E)}$ of the fields also shows significant dependence on the local pixel morphology, where the value of the phase difference between two kinds of pixels is $\mathrm{\pi}$, confirming the prediction by the scalar-based theory (see the texts below) [29]. Compared to Figs. 1(a) and 1(a’), Figs. 1(b) and 1(b’) clearly show that both the intensity and the phase of the transmission fields exhibit across-pixel continuity and display the sub-pixel profile for the DOE with $\mathrm{\Lambda_{ph}/\lambda=0.6}$. The assumption of uniformity and local-pixel-dependence of optical fields made by the scalar-based theory, i.e., $\mathrm{|E|=1}$ and $\mathrm{arg(E)=2\pi n_{1}h/\lambda}$ for the transparent pixel, and $\mathrm{|E|=1}$ and $\mathrm{arg(E)=2\pi n_{2}h/\lambda}$ for the pink pixel, is hence considerably invalid. Such results may also imply that these single-value representation of optical fields per pixel by the scalar-based theory can be problematic in the regime $\mathrm{\Lambda_{ph}/\lambda\leq 1}$.

To explore the generality of vectorial effects for systems studied, we statistically analyze fields of a specified (central) pixel inside DOE as a function of pixel size, by means of an ensemble of 20000 specimens per $\mathrm{\Lambda_{ph}}$ value. Similar to the schematic DOE in Fig. 1(c), each specimen has its plane dimension $3.5\mathrm{\lambda\times 3.5\lambda}$ and is assigned with stochastic pixel morphology except for the specified (central) one. Here, the FDTD method is adopted for giving convincing conclusions. Numerical results for amplitudes of different field components on the exiting surface (at z=$\mathrm{\lambda+\delta_{h}}$) are plotted in Fig. 2(a), provided that the incident field has x polarization. The curve of phase of the field’s x component is shown in Fig. 2(b), where the curves for the other two components are not displayed, owing to their simple uniformly random distribution over the 0-2$\pi$ range. In principle, Fig. 2 portrays the properties of optical fields from the regime $\mathrm{\Lambda_{ph}/\lambda\gg 1}$ governed by the scalar-based theory (ST), across the regime $0\ll\mathrm{\Lambda_{ph}/\lambda\leq 1}$ belonging to vectorial optics, and to the regime $\mathrm{\Lambda_{ph}/\lambda\ll 1}$ conforming to the effective medium theory (EMT)[16]. As expected, numerical results show that the amplitude and phase of fields exhibit asymptotic convergence toward constants $\mathrm{E_{x,EMT}}$ and $\phi_{\mathrm{E_{x,EMT}}}$ by EMT at $\mathrm{\Lambda_{ph}/\lambda\ll 1}$, as well as toward constants $\mathrm{E_{x,ST}}$ and $\mathrm{\phi_{E_{x,ST}}}$ by ST at $\mathrm{\Lambda_{ph}/\lambda\gg 1}$, respectively. Considerable depolarization and dephasing of fields, however, occur in the in-between regime and hence demand vector-based theories for different requirements.

2.2 K-series approximation method to compute pixel fields $\mathrm{\mathbf{P}(x,y)}$ by the vectorial model

To include vectorial effects into designs of finite-size (mm- to cm-scale) devices (with structure parameter $\mathrm{\Lambda_{ph}/\lambda\leq 1}$) for practical applications, the two most popular rigorous electromagnetic theories, FDTD and RCWA, have been adopted for various attempts as discussed in Sec. 1. Unfortunately, the FDTD method signifies highly-challenging issues, since its computation time scales are the fourth order of the simulation (spatial) domain size N (e.g., $\mathrm{(10N)^{4}\approx 10^{(1+5)\times 4}}$ at cm-scale), and the memory requirements are the third order [30]. The RCWA method has the same issues due to its computation time, as well as memory requirements, being fourth order of the simulation (momentum space) domain size K [31] (e.g., $\mathrm{K^{4}\approx 10^{5\times 4}}$ at cm-scale). To efficiently design optics of DOEs, on the basis of the semi-analytical RCWA framework, this subsection presents a k-series approximation method associated with optional parameters, the rank of neighboring interferences r, and the sampling rate s, for selecting desired degrees of accuracy at the cost of computational resources. In principle, the method scales the computation time and memory requirement down to the level of $\mathrm{N^{2}(2r+1)^{2}s^{4}}$; e.g., about $\mathrm{10^{11}}$ at condition r=s=1 at cm-scale, in comparison with the level of $\mathrm{N^{2}}\simeq 10^{10}$ by the scalar-based method [29]. Here, the translation symmetry has been considered to estimate the computational cost for KSAM. We note that, implemented with lookup table techniques, the proposed method can further scale the computation time down to the level of $\mathrm{N^{2}}$; whereas too large values of r and s can burden this method with high computational costs beyond FDTD and RCWA. The underlying criteria are described in detail in the following paragraphs.

For a binary DOE (one stacking layer), its permittivity and the incident field can be represented in a Fourier series of spatial harmonics [23] according to

[TABLE]

in which the wave vector components $\mathrm{k_{xu}}$ and $\mathrm{k_{yv}}$ arise from phase matching and the Floquet conditions and are given by

[TABLE]

Here, $\mathrm{\Lambda}_{x}$ and $\mathrm{\Lambda}_{y}$ are the periodicity along the x and y directions, respectively; $\mathrm{\lambda}$ is the vacuum wavelength of the optical wave; $\mathrm{k_{0}=2\pi/\lambda}$ is the free space wave number; $\mathrm{n}$ is the refractive index of the incident/outgoing medium ($\mathrm{n}=1$ for the air in this work); $\mathrm{\theta}$ is the polar angle; and $\mathrm{\varphi}$ is the azimuth angle of the incident k-vector. We notice that the vector variables are in bold for convenience. On the basis of the spatial harmonics, Maxwell’s equations can be expressed as an infinite set of first-order differential equations:

[TABLE]

Together with specific incidence conditions, Eqs. (8)-(11) can directly solve the vectorial fields $\mathrm{\mathbf{e}}$ and $\mathrm{\mathbf{h}}$ by means of eigenproblem algebra and can be easily extended to $\mathrm{\ell}$-level ($\mathrm{\ell}=2$ for the binary one) DOEs [29, 23]. For computationally efficient designs of DOEs, this work computes vectorial fields of an individual pixel under the pre-condition that the (singlepixel) field includes accurate characteristics of across-pixel interferences and sub-pixel fringes. In this way, KSAM can appeal to very finite k-series for still predictive computations, since the interference fields occur only among finite nearby pixels and the sub-pixel fringes are moderate under the condition $\mathrm{0\ll\Lambda_{ph}/\lambda\leq 1}$. By the reciprocity relations in Fourier algebra, we correlate the order of k-series with spatial parameters- i.e., the rank of neighboring interferences r, the sampling rate s, and the DOE depth h- for easily selecting the desired degrees of accuracy.

Assuming negligible multiple reflections of electromagnetic waves, this work introduces an effective dimension $\mathrm{\Lambda_{eff}\equiv(2r+1)\Lambda_{ph}}$ to compute the fields for an assigned pixel, with the condition that the set-up can reproduce essential interferences’ effects and simultaneously eliminate disturbances from non-physical periodicity in RCWA algebra. Moreover, since the spatial frequency $\mathrm{\varepsilon_{u-m,v-n}}$ of every stacking layer in DOEs can iteratively transfer wave vectors $\mathrm{\mathbf{e}_{mn}(\mathbf{h}_{mn})}$ into different diffraction components $\mathrm{\mathbf{e}_{uv}(\mathbf{h}_{uv})}$, as indicated in Eqs. (8)-(11), simple criteria for spatial parameters can be given by $\mathrm{r\Lambda_{ph}\geq hk_{\parallel,R}/k_{z,R}}$ and $\mathrm{s\geq\ell-1}$. Here $\mathrm{k_{z,R}}$ is the selected highest-order k vector with a real-number z component, while $\mathrm{k_{\parallel,R}=(k_{x,R}^{2}+k_{y,R}^{2})^{1/2}}$ is the in-plane component of the k vector defined for $\mathrm{k_{z,R}}$. The total components of k-series are selected in the range $\mathrm{-(2r+1)s/2\leq u,v\leq(2r+1)s/2-1}$ with Eqs. (5)-(2.2) at $\mathrm{\Lambda_{x(y)}=\Lambda_{eff}}$ in this work. The graphical interpretations of the optional parameters r and s are illustrated in Fig. 3.

In Fig. 3(a) the rank of neighboring interferences r defines the maximal pixel range of pixel-pixel interferences by $(\Lambda_{\mathrm{eff}}/\Lambda_{\mathrm{ph}}-1)/2$, where the size of an equivalent-isolated pixel $\mathrm{\Lambda_{eff}}$ should be infinite in theory for accuracy, but is set as a finite value by assuming insignificant interactions of fields among distant pixels. The computation hence requires only a finite morphology of nearby pixels within $\Lambda_{\mathrm{eff}}$ to determine the local field in $\Lambda_{\mathrm{ph}}$, but still can keep predictive calculations. In Fig. 3(b) the sampling rate s defines the number of sampling points per pixel in one dimension. In this way, the k-series can represent an accurate sub-pixel profile of fields and enable occurrences of higher-order diffractions by reciprocity relations. With an adequate set-up of r and s values at given h, this work fulfills more efficient and accurate vectorial modeling of near fields in DOEs, in place of the typical wave front function in the scalar-based theory. Figure 4 illustrates schematic examples for computing the fields of the (red) pixel P(0,0) at different conditions (r,s) for KSAM, in which the lower part indicates the pixel morphology $\mathrm{\Upsilon}$ (with total dimension $\mathrm{\Lambda_{eff}}$ and $\mathrm{\Lambda_{ph}/\lambda=1}$), and the upper part shows the profiles of fields $\mathrm{|E_{x}(x,y)|}$ at $\mathrm{z=\lambda+\delta_{h}}$. Numerical results state that one can perform predictive optical designs with proper dimension r (comparing Figs. 4(a) and 4(b)) and sampling data s (comparing Figs. 4(b) and 4(c)), at an affordable computational cost ($\mathrm{92\%}$ reduction in computational time by a comparison of Figs. 4(a) and 4(c)). Here, the spatial domain size N is 1 for determining the (red) single-pixel fields by KSAM.

To investigate the validity for general cases, we re-evaluate the analyses in Fig. 2 by KSAM and compare results with those by FDTD. Figure 5(a) investigates the average amplitude of different field components at the specified (central) pixel using r=s=2. Figure 5(b) plots the curve for the corresponding phase of x components of fields. Curves by FDTD are drawn by dotted lines for comparison. The asymptotic constant values based on the EMT and ST methods are denoted as horizontal dotted lines near the corresponding regimes, respectively. Numerical results in Fig. 5 display qualitatively similar variations of fields to that by FDTD when varying $\Lambda_{\mathrm{ph}}$. With higher-order (r=s=4) approximations, Fig. 6 shows quantitatively similar curves to that by FDTD and clearly predict more accurate vectorial effects (e.g., depolarization and dephasing of fields) for systems studied.

2.3 Vector-based simulated annealing algorithm for designing diffractive optical elements

To next interpret the process for optical designs, this work exemplifies a Fourier-transformed DOE, where the coordinates of the image plane are denoted ($\mathrm{x_{I}}$,$\mathrm{y_{I}}$) and those of the exiting plane on DOE are denoted ($\mathrm{x_{D,s}}$,$\mathrm{y_{D,s}}$). DOE is assumed to be illuminated by a plane wave. The reconstructed image is derived in the focal (image) plane of the Fourier lens.

During the design process, trial morphologies $\mathrm{\Upsilon^{\prime}(x_{D},y_{D})}$ are specifically first generated with a simulated annealing algorithm, and the pixel-field profiles $\mathrm{\mathbf{P}^{\prime}(\mathrm{x_{D,s}},\mathrm{y_{D,s}})}$ are then computed by means of KSAM. Here, the subscript s represents the index of sampling points per pixel. The DOE quality hence is scored by a cost function associated with the diffraction image of $\mathrm{\mathbf{P}^{\prime}}$ and the target image, in which the diffraction fields $\mathrm{R_{x(y)}}$ in the image plane are given by

[TABLE]

Here, $\mathrm{F^{-1}}$ is the inverse Fourier transform, and the $\mathrm{R_{z}}$ component is obtained with the relation $\mathrm{\mathbf{k}(x_{I},y_{I})\cdot\mathbf{R}(x_{I},y_{I})=0}$. Moreover, DOE is considered as an optical element of $\mathrm{N^{2}}$ square pixels, and the pixel-field profile $\mathrm{\mathbf{P}}$ is an $\mathrm{(N\times s)^{2}}$ array function. In this work $\mathrm{x_{I}}$ and $\mathrm{y_{I}}$ are simply assigned as the lowest $\mathrm{N^{2}}$-order wave vectors of the total $\mathrm{(N\times s)^{2}}$ ones, responding to an $\mathrm{N^{2}}$-pixel target image $\mathrm{I_{0}}$. The intensity I of the reconstructed image is given by:

[TABLE]

The cost or energy function to score performance of reconstruction is defined as the mean-square error between the reconstructed image and the target images $\mathrm{I_{0}}$, presented by:

[TABLE]

Here, $\alpha$ is a scale factor given by:

[TABLE]

We note that the cost function is obviously zero when the reconstructed image equals the desired image.

The vector-based simulated annealing algorithm can be summarized by the following steps, and the flow chart is plotted in Fig. 7.

The pixel morphology $\mathrm{\Upsilon(x_{D},y_{D})}$ of DOE is initialized; e.g., according to the discretization of phase distribution of the Fourier transform of a desired image $\mathrm{I_{0}}$. The corresponding pixel-field profile $\mathrm{\mathbf{P}(x_{D,s},y_{D,s})}$ is then initially built by KSAM. The temperature parameter $\mathrm{T_{0}}$ of the simulated annealing process is initially set at a relatively high value, enabling a large range of perturbation probability. 2. 2.

Each pixel $\mathrm{\Upsilon(x_{Dr},y_{Dr})}$ of DOE is successively selected for fluctuation during the pixel cycle ($\mathrm{x_{Dr}=1:N}$ and $\mathrm{y_{Dr}=1:N}$). The morphology (or depth) of this selected pixel is reset by a random (trial) candidate of depth levels. The fields of pixel $\mathrm{\mathbf{P}^{\prime}(i,j)}$ in the range ($\mathrm{x_{D}-r\leq i\leq x_{D}+r}$,$\mathrm{y_{D}-r\leq j\leq y_{D}+r}$) are then updated by KSAM, where these pixels belonging to device boundaries are assumed to follow a periodic arrangement for simplification. We note that updating multiple pixel fields $\mathrm{\mathbf{P}^{\prime}(i,j)}$ for a single trial pixel $\mathrm{\Upsilon(x_{Dr},y_{Dr})}$ here can bring forth a factor of $\mathrm{(2r+1)^{2}s^{2}}$ for the computation time. 3. 3.

Calculate the diffraction fields by Eq. (2.3). The difference in the cost function $\mathrm{\Delta m}$=$\mathrm{m_{new}}$-$\mathrm{m_{now}}$ is then calculated by Eqs. (13)-(15). If $\mathrm{\Delta m}$ is negative, then the new pixel morphology is accepted; otherwise, acceptance or rejection is determined by the Boltzmann probability:

[TABLE]

where a random number $\mathrm{0\leq\zeta\leq 1}$ is first generated, and the new (trial) pixel is accepted if $\mathrm{\zeta<\Theta(\Delta m,T)}$ occurs and is rejected otherwise. 4. 4.

Repeat steps 2 and 3 for the pixel cycle ($\mathrm{x_{Dr}}$ and $\mathrm{y_{Dr}}$) and the annealing cycle ($\mathrm{n=1:N_{A}}$); except when the criterion of convergence is achieved, where T in Eq. (16) is made lower with the number of annealing cycles by:

[TABLE]

An alternate Fourier algorithm for free-space propagation of vectorial optical fields, compared to Eq. (2.3), can be seen in [32]. The Matlab code concerning the first-order algorithm (r=1 and s=1) for designs of square DOEs is available in MATLAB Central File Exchange [33].

2.4 Two-stage optimization process for high-quality DOEs

To achieve high-quality DOEs, we introduce complementary optimization processes for deviations by hypothetical approximations. Specifically, for systems with microstructures at the subwavelength scale, field errors are primarily attributed to two issues: vectorial effects and approximation by the truncation of k-series. Thus, a two-stage optimization process is proposed in this work to deal with these two issues. To enable fast estimations, we employ a scaled-down sample of the desired product as a temporary prototype during the two-stage process. A schematic flow diagram is plotted in Fig. 8. Referring to Fig. 2, the vectorial optical fields in the studied regime $\mathrm{0\ll\Lambda_{ph}/\lambda\leq 1}$ noticeably depart from the prediction by the scalar-based theory. The optimized DOE depth, binary DOE for instance, could hence differ from the ideal value as well. For this reason, it is necessary to identify the optimal depth h by the vector-based theory in the first stage. By applying separate DOE designs $\mathrm{\Upsilon^{\prime}}$ at different depth value, the highest-quality one $\mathrm{\Upsilon^{\prime}_{d}}$ and the optimized design depth $\mathrm{h_{d}}$ are determined, as processes (2) and (3′) in Fig. 8. Here, the popular specification by the 0-order diffraction is adopted as the quality criterion. In the second stage, since the transmission fields by low-order KSAM deviate from that by the rigorous FDTD method (see Figs. 5 and 6), the design depth $\mathrm{h_{d}}$ in theory can be inconsistent with the optimized value $\mathrm{h_{p}}$ for products. To determine the optimal $\mathrm{h_{p}}$ value, the design $\mathrm{\Upsilon^{\prime}_{d}}$ is scored as a function of depth h by FDTD. It is concluded that the final DOE product $\mathrm{\Upsilon_{p}}$ is designed with v-SAA at $\mathrm{h_{d}}$ and is manufactured with the set-up of depth $\mathrm{h_{p}}$ for optimal performances. One can take it for granted that the value $\mathrm{h_{d}}$ will asymptotically converge to the value $\mathrm{h_{p}}$ in the high-order limit of the k-series approximation.

3 Design results and analysis

To evaluate the feasibility of the proposed algorithm, a beam shaper is adopted as an exemplary target. The simulation conditions are set by the plane-wave incidence, the wavelength $\mathrm{\lambda=1.0}\ \mathrm{\mu m}$, the incident angle $\mathrm{\theta=0^{\circ}}$, TM polarization, the refractive index of air $\mathrm{n_{1}=1.0}$, the refractive index of DOE material $\mathrm{n_{2}=1.5}$, the physical pixel size $\mathrm{\Lambda_{ph}=0.6\mu m}$, the rank of neighboring interferences r=1, and the sampling rate s=1 position per pixel, for designing a $\mathrm{201\times 201}$ pixel array of binary DOEs. The default depth is set as $\mathrm{h_{def}=1.0\lambda}$ by s-SAA, and its optimal value is re-evaluated by the two-stage optimization process as in Fig. 9. A 3D schematic of DOEs is shown in Fig. 1(c). The full structure of DOEs is represented by an image matrix, as in Figs. 10(a) and 11(a), where the cyan area indicates the material medium ($\mathrm{n_{2}}$), and the blue area is the air-filling region ($\mathrm{n_{1}}$). The beam shaper target is indicated as the inset mesh in Figs. 10(c) and 11(c).

Before generating the full $\mathrm{201\times 201}$ design, a temporary $\mathrm{101\times 101}$ DOE is employed for speeding up the two-stage optimization process, so as to quickly determine the optimal design depth $\mathrm{h_{d}}$ and product depth $\mathrm{h_{p}}$. At the first stage, with varying depth parameter h, a different DOE morphology is separately generated, and its intensity of zero-order diffraction is recorded for a quality score. Numerical results in Fig. 9(a) show that the optimal design depths are $\mathrm{h_{d,sSAA}=\lambda}$ and $\mathrm{h_{d,vSAA}=1.05\lambda}$, corresponding to the binary morphologies $\mathrm{\Upsilon^{\prime}_{d,sSAA}}$ and $\mathrm{\Upsilon^{\prime}_{d,vSAA}}$ by s-SAA and v-SAA models, respectively. At the second stage, with given designs ($\mathrm{\Upsilon^{\prime}_{d,sSAA}}$ and $\mathrm{\Upsilon^{\prime}_{d,vSAA}}$), the zero-order intensities at different set-up of depths are directly estimated by FDTD, in order to determine the optimal depth $\mathrm{h_{p}}$ for products. Figure 9(b) shows that the optimal product depths are $\mathrm{h_{p,sSAA}=h_{p,vSAA}=1.15\lambda}$ for this case. Here, the large zero-order intensity signifies the low-quality DOEs that are designed by the s-SAA model at a subwavelength feature scale.

With the known optimal depth parameter $\mathrm{h_{d}}$, the full $\mathrm{201\times 201}$ morphology of DOE can then be generated by SAA. Figure 10(a) illustrates the binary morphology of DOE designed by the s-SAA model with $\mathrm{h_{d}=1.0\lambda}$. The 4-thread computation over 100 annealing cycles runs in 160.9 seconds, which converges the variation of the cost function to below $\mathrm{0.01\%}$ per iteration. The $\mathrm{201\times 201}$ morphology of DOE designed by the v-SAA(r=1,s=1) model with $\mathrm{h_{d}=1.05\lambda}$ is depicted in Fig. 11(a). Similarly, the 4-thread optimization computation over 100 cycles runs in 1336.9 seconds, plus an extra time of 1.4 seconds for building the lookup table having $\mathrm{\ell^{(2r+1)^{2}}s^{2}}$-size. To validate the designs, with given $\mathrm{h_{p}=1.15\lambda}$, the reconstructed images of DOEs are directly simulated by means of well-known softwares [34, 35] that use the rigorous electromagnetic theory FDTD. Figure 10(b) shows the reconstructed image that corresponds to the design in Fig. 10(a), and its height plot in lower resolution is shown in Fig. 10(b’). The target pattern is put as an inset mesh in Fig. 10(b’) for comparison. Here, rather than the typical intensity image $\mathrm{I=|E|^{2}}$, the reconstructed image is presented by the amplitude value of fields $\mathrm{|E|}$ for clarity. The FDTD program by 4-threads estimates the optics of DOE in 23389.1 seconds, in comparison with 1.6 seconds by s-SAA and 13.4 seconds by v-SAA. Numerical results indicate that the DOE designed by the s-SAA model exhibits low diffraction efficiencies, owing to the overwhelming zero-order intensity. Numerical results for the design by v-SAA are reported in Fig. 11. The corresponding reconstructed images are illustrated in Fig. 11(b) and 11(b’). Accordingly, we find that diffraction efficiencies as well as the design quality can be significantly enhanced by the vector-based method. Beyond the mentioned cases, designs with a finite device extent of, multi-level ($\mathrm{\ell\geq 3}$), and/or multiple functional diffractive optical elements could also be generated by the presented model.

4 Conclusion

This research has described in detail a vector-based optimization tool for designing DOEs with subwavelength feature sizes. By introducing an adequate selection scheme of k-series, the algorithm can perform computationally efficient yet with predictive computations. Together with an optimization of the geometrical degrees of freedom (i.e., DOE depth h) as compensation for errors from the truncation of k-series, this work presents DOEs with diffraction efficiencies close to the theoretical limit. Furthermore, after implementing with lookup table techniques to pre-calculate pixel fields under ergodic couplings with nearby ones, the vector-based algorithm can achieve large-extent (mm- to cm-scale) designs in time frames comparable to those of scalar-based models.

\bmsection

Acknowledgments This work was supported by Institute of Nuclear Energy Research (INER) under project AIE010403.

\bmsection

Disclosures The authors declare no conflicts of interest.

\bmsection

Data availability The Matlab code concerning the first-order algorithm (r=1 and s=1) in this paper is available in Ref. [33].