Inverse design of diffractive optical elements using step-transition perturbation approach

Diffractive optical elements are ultra-thin optical components required for a variety of applications because of their high design flexibility. We introduce a gradient-based optimization method based on a steptransition perturbation approach which is an efficient approximation method using local field perturbations due to sharp surface profile transitions. Step-transition perturbation approach be available to calculate the gradient of figure of merit straightforwardly, we implemented optimization method based on this gradient. This fast and accurate inverse design creates binary (2-level) diffractive elements with small features generating the wide angle beam arrays. The results of the experimental characterization confirm that the optimization based on the perturbationmethod is valid for 1-to-117 fan-out grating generating beam pattern of linear array.


Introduction
Diffractive optical elements (DOEs) are used for a variety of optical systems because of their compact size, high design flexibility, and ease of mass production. A good example of this type of device is fan-out grating, often also referred to as diffractive beam-splitter which creates multi spots by deflecting an incident light into different diffraction orders.
Applications of fan-out elements include multi-focal microscopy [1,2], camera calibration [3], optical system distortion measurement [4], optical interconnects [5], and structured light projectors [6,7]. In particular, wide angle DOEs are used in promising field of applications [7][8][9] with recent progress in fabrication technology realizing nanoscale features. When only small diffraction angles are required, the iterative Fourier transform algorithm (IFTA) [10][11][12] based on the thin element approximation (TEA) [13] is widely and successfully used for the design of the microstructure surface of the DOEs. However, this approach suffers from several severe shortcomings and is no more valid when larger angles are required [14]. One of the most severe problems is the insufficient modeling of the light field transmission through the DOE by the TEA [15]. A precise modeling can be obtained by rigorous electromagnetic diffraction theory, such as rigorous coupledwave analysis (RCWA) [16][17][18]. However, the use of parametric optimization based on the rigorous analysis is often computationally heavy because the gratings with many parameters are leading to high-dimensional optimization problems. To overcome this limitation, we employ gradient-based algorithms [19] based on the step-transition perturbation approach (STPA) [20]. In gradient-based optimization of optical elements, within the overall design parameter space, one iteratively searches the maximum or minimum point based on calculated gradient [21,22]. STPA enables rapid approximate calculation of diffraction intensity of wide angle DOEs, which is an approximate method based on local field perturbations generated by sharp transitions of the surface profile of diffractive elements. We can describe analytically the gradient of diffraction efficiency with respect to the design variables of optical elements using STPA. It thus is feasible to calculate the gradient straightforwardly with accuracy as much as the approach based on the rigorous electromagnetic diffraction theory if most of the features of the structure are bigger than the wavelength of the incident light.
Here, we introduce a design approach based on STPA for wide angle DOEs, yielding improvements in the uniformity of the created patterns while maintaining the total diffraction efficiency. We have focused our efforts on designing binary (i.e., 2-level) micro-structures because they are most easy to fabricate and thus obviously are very attractive for optical systems.
A schematic of DOE generating linear spot array is shown in Figure 1(a). In order to obtain a starting condition for our optimization, a one-dimensional (1D) diffractive fan-out DOE is designed by TEA-based IFTA (see Figure 1(b) inset). Figure 1(b) present the performance of this initial grating as a function of diffraction angle. Obviously, the performance of this initial fan-out element is very unsatisfactory with respect to uniformity even when the maximum diffraction angle is over 7°. We then apply our method to optimize large-angle 1D diffractive beam splitters but also compare to the results from parametric optimization based on RCWA. We measured the performance of fabricated samples based on the optimized designs and compared to calculated diffraction efficiency using RCWA simulation. Based on this analysis, we can show the excellent strengths of our design method.

Inverse design methods
An important aspect of the optimization process is the parametrization used to describe the shape of the optical elements, which can significantly affect the performance and computational cost. Figure 1 illustrates an example of a 1D binary phase grating profile with 2K transitions in position x k within a single grating period. We use these positions of transition points as the set of design parameters x = [x 1 ,…,x k ,…,x 2k ] and define the figure of merit (FOM) to optimize DOEs creating diffraction pattern with uniform intensity distribution: where F represents the difference between the calculated diffraction efficiency η m and the target diffraction efficiency η obj in diffraction orders. The gradient of the FOM with respect to transition positions ∇ x F is crucial in determining the search direction to optima. For example, if the total number of transitions 2K is large, it may easily become computationally heavy to calculate the gradient by RCWA analysis. The STPA, however, allows expressing the variation for a diffraction efficiency with respect to transition positions as an analytical solution so that it can calculate the gradient straightforwardly.
It has been reported by T. Vallius et al. [20] that in fact the approximated method based on local field perturbations from sharp step-transitions enables rapid calculation of diffraction patterns of DOEs in the nonparaxial domain. The perturbations are observed in the field distribution directly after sharp vertical transitions of binary gratings. The TEA calculation, however, yields a constant amplitude and phase. This omission of perturbations in TEA makes computing inaccurate gratings with wavelength-scale structures, i.e., the gratings creating the wide angle arrays. Thus, we can accurately calculate the diffraction efficiency using the model which combines the TEA with field disturbances caused by sharp transitions in the surface profile calculated by RCWA. We define the field perturbation behind the kth sharp transition located at the point x k in the surface profile as where U R k (x) and U T k (x) are field calculated by RCWA and TEA, respectively and Δ T is the truncation parameter that is chosen 10λ in the calculations [20]. In Figure 2(a), (b), the amplitude and phase of the field distribution directly after an isolated step transition determined by TEA and RCWA is presented. The field perturbations of binary gratings consist of only two kinds of oscillation corresponding to left-side and right-side transition point in a ridge. Therefore, the constructed field behind binary grating with many transition points is described by the x-axis shifts of the two field perturbations p 1 (x) and p 2 (x) in the following expression: where 2K is the total number of the transitions. The amplitude and phase of the field perturbation of right-side of a ringe p 2 (x) is represented in Figure 2(c), (d). The diffraction amplitude of mth order in far field is given by where Λ is the grating period and T m and D m is the Fourier coefficient of the field calculated by TEA and a field perturbation contribution, respectively.
where the Fourier coefficient P m of field perturbation p 1 (x) is expressed as The Fourier coefficient of p 2 (x) is P −m in Eq. (6) because the p 2 (x) is an even function of p 1 (x). We found D m in Eq. (6) using the Fourier shifting theorem [23]. Up to here, T. Vallius et al. claimed that this is an efficient computation method compared with the calculation of RCWA to the entire profile. Once the Fourier coefficient of the steptransition perturbation P m and P −m is calculated and no further RCWA calculations are necessary.
We furthermore focus on the fact that this Fourierdomain contribution from step transition P m don't contains explicit dependence on transition point x k . This point is highly useful when calculating the gradient of diffraction efficiencies with respect to transitions positions ] to optimize the structures. To find these derivatives, we apply chain rule when differentiating the FOM F (x): where the first term ∂F ∂η m is easily calculated by using Eq. (1) and the second term is also expressed by an analytical equation because because P m and P −m don't include the dependence on the position of transition point x k . Additional details on this are described in the Appendix. It is feasible to calculate the gradient straightforwardly with accuracy as much as the approach based on the rigorous method if most of the features of the structure are bigger than the wavelength of the incident light. The obtained gradient was used in optimization based on the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [24,25].

Simulation results
Using the proposed optimization approach, we can design various multi-spot array generators. In general, diffractive beam splitter creating larger number of spots require more complex structure, i.e. gratings with many features. To verify our method is valid in high dimensional optimization problems, we show the optimization results of fan-out grating generating many spots, for instance, 1-to-117 diffractive beam splitter.
To evaluate DOEs, we use two different metrics which are uniformity error (UE) and normalized root-meansquare error (NRMS) σ.
where η max and η min represent the maximal and minimum diffraction intensity and η m is diffraction efficiency in orders and η obj is target diffraction efficiency and M is the total number of diffraction orders. Lower values of both UE and NRMS indicate less residual variance so that our objective is to minimize UE and NRMS of a DOE design given uniform diffraction efficiency distribution.
To apply the optimization method, we prepared surface profiles of 1D fan-out grating designed by IFTA. The fused silica (SiO 2 ) was selected as the material. The refractive index of SiO 2 is assumed as n 2 = 1.46. Transverse electric (TE)-polarized (i.e., E-field component along the yaxis) monochromatic light with a wavelength of λ = 633 nm is an incident plane wave from the substrate side with normal incidence angle. The depth of the grating was selected as d = 692 nm and the grating period is 200 µm. Thus, the maximal diffraction angle of 1-to-117 diffractive beam splitter are about 11°at 58th order from 0th order.
To optimize these 1-to-117 diffractive beam splitters, we use our figure of merit as in Eq. (1) with the uniform intensity distribution of target efficiency η obj and find the local optima using the L-BFGS algorithm with the gradient calculated based on STPA. The uniformity of beam array created by elements designed from IFTA based on TEA, followed by optimization, are plotted in Figure 3 with different metrics. We also plot together with the uniformity of final design after gradient-based optimization based on RCWA, in this case, the gradient calculated by brute-force approaches. In other words, normally around 60 times (i.e., the number of transition points 2K + 1) RCWA simulation is required in an iteration during the optimization. We compared the optimized results by gradient-based on STPA and RCWA. For an accurate comparison, all diffraction efficiencies of final designs are calculated by RCWA. In most cases, the uniformity of these final elements is significantly improved and the uniformity of final design optimized based on STPA are as good as those of optimized based on RCWA. However, the performance of optimization based on STPA is much better than based on RCWA in terms of computation effort. The simulation and optimization steps were written using MATLAB scripts, and the optimization process took less than 20 s using gradientbased optimization by STPA, while taking over 6 h using gradient-based optimization by RCWA on a machine with 3.60 GHz clock rate and 32 GB RAM. During the optimization, the diffraction pattern for calculating UE and NRMS were evaluated with RCWA solver RETICOLO [26].
To observe the changes of the FOM and transition positions during the optimization, we represent one optimized 1-to-117 diffractive beam splitter in Figure 4.   Figure 4(b). The total number of transitions is 66 and the average change of transition points is around 300 nm after optimization. The simulated diffraction efficiency distributions of DOEs after optimization is shown in Figure 4(c). We calculated the total diffraction efficiency, UE, and NRMS of optimized diffractive beam splitters. The total diffraction efficiency of 117 spots of optimized DOE is 77.35% and UE from 38.68 to 10.79% and NRMS from 12.16 to 04.18%, through gradient-based optimization using STPA. The surface profile of optimized design which has critical dimension (CD) (i.e., minimum feature size) is 700 nm and fill factor is 51.16% is represented in Figure 4(c) inset.

Experimental results
The diffractive beam splitters were fabricated by direct laser writing to obtain SiO 2 binary surface relief structures. The elements are optically characterized using a  TE-polarized 636 nm wavelength beam from a diode laser. We detect the diffracted light beams using a mobile singlepixel detector with a high dynamic range. In Figure 5, a detector with a pinhole aperture is mounted on a translation stage under computer control. By scanning the detector to the center of each of the spots, it is possible to measure the power contained in each of the spots, i.e. diffraction orders in the output array.
To focus on both the simulation and experiment to facilitate a quantitative comparison, we applied loss caused by Fresnel reflection from the interface between air and SiO 2 substrate to simulate the overall efficiency of DOEs. The comparison between theoretical and experimental diffraction efficiencies are presented in Figure 6. We represent the total diffraction efficiency, UE, and NRMS of simulated and measured one in Table 1.
The experimental data show that the DOE operates with high-performance. The UE and NRMS of beam splitters are 21.42 and 8.07%, respectively. For an accurate comparison between theoretical and measured results, we analyze the correlation of these data using mean absolute percentage deviation (MAPD) as a ratio defined by the formula:   where η S m , η E m are simulated and experimental efficiency in (m)th diffraction orders and M is the total number of diffraction orders. The MAPD of 1-to-117 beam splitters are calculated to the 8.15%, which shows excellent reproducibility of the simulated results in a quantitative manner. The only noticeable deviation in the measurement is a small mismatch of diffraction efficiency in a few orders due to minor fabrication errors. In general the diffraction efficiency in orders often strongly depends on the errors in fabrication processes, e.g., etching depth, feature width, slope steepness, and feature rounding. Nevertheless, the fabricated samples based on optimized design overall display experimental performances which are better than the theoretical performances of initial designs before optimization.

Conclusion
In summary, we utilized the STPA in optimizing the optical elements, which is able to create wide angle diffractive optical elements at a very low computational cost. We explored properties of the optimization method, such as efficient computation for the gradient of the target function with respect to transition positions with Fourier-domain local field perturbation. As a case study, we applied gradient-based optimization with STPA to 1-117 beam splitter with a non-paraxial diffraction angle, i.e., maximal diffraction angle is 11°from the center, respectively. The optimized beam splitter show a considerable improvement of uniformity while maintaining the initial diffraction efficiency. The experimental results obtained by the illumination of the fabricated optical elements using a laser of 635 nm wavelength with a normal incidence have been compared with the numerical results. Numerical simulation and experimental results were found to be in good agreement and our optimization method can be considered proven to be an effective design tool for wide angle diffractive beam splitters.
Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.
Research funding: Horizon 2020 Framework Programme (675745). Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

Appendix
In order to calculate the gradient of figure of merit in Eq. (8), we calculated the derivatives of the diffraction efficiencies where the diffraction efficiency η m is a function with respect to transition point x k in binary grating. Thus we can partially differentiate each term of η m with respect to x k . when m ≠ 0, we can express the derivatives as where Φ m = sin(Δϕ/2)/πm (14a) The Δϕ is the difference between phase ϕ 1 and ϕ 2 which are the phase of an electric field in the air and dielectric material, respectively and 2K is the number of transition point in structure. The Fourier coefficients P m and P −m of field perturbation are given by Eq.