1. Introduction

Topology optimization is the most powerful tool for structural optimization because it extends the design space beyond that of size and shape optimization, and it can also provide useful designs for high-performance structures that implement new targets. More than 30 years have passed since Bendsøe and Kikuchi [1] first proposed the concept of topology optimization in the field of computational mechanics. During this period, many pioneering studies were published that promoted the development of related theories and technologies, and spurred the use of practical topology optimization methods, such as the density-based method [2,3], the level set method [4,5], and so on [6,7].

As topology optimization gradually matures, scholars have begun to extend its applications from computational mechanics to more complex problems. Some innovative work on phononic crystals [8,9], photonic crystals [10,11], multiscale systems [12,13], acoustic vibration [14], and metalens optimization [15,16] continue to be proposed.

Recently, the inverse design computational method has received significant attention in nanophotonic devices, and a series of high-efficiency designs have been established based on the derivation of adjoint variables [17,18]. Lin et al. [19] proposed a general inverse design framework for metasurfaces that can automatically discover highly complex multi-layer superstructures by using density-based topology optimization. Piggott et al. [20] designed a silicon wavelength demultiplexer and successfully manufactured it to demonstrate its effectiveness. Chung et al. [21] applied inverse design to broadband and achromatic metalenses, showing the highest theoretical efficiencies in low numerical apertures. Lin et al. [22] used topology optimization to solve an optical metasurface optimization problem with 10⁵∼10⁶ degrees of freedom and obtained a novel design. Some educational articles have also emerged on the optimization of optical topologies, such as online databases [23], code writing [24], and software applications [15].

In the design of nanophotonic devices, three-dimensional (3D) laser nanolithography is among several technologies with the highest degree of manufacturing freedom. However, this processing method is still affected by the inherent shortcomings of additive manufacturing technology, that is, the need for supporting structures. This requires additional processing constraints [25] or length-scale control [26]. Augenstein and Rockstuhl [27] proposed design schemes that combine compliance and photonic devices performance from multi-objective design. Wang et al. [28] combined topology optimization and length-scale control methods to design manufacturable photonic crystal cavities.

The recently proposed material-field series-expansion (MFSE) method [29,30], which inherently avoids the checkerboard patterns and mesh dependency by adopting the spatially dependency material-field definition, can be applied to the optical metalens problem to achieve some competitive topology designs. Because the optical metalens is a typical nanophotonic device that should have controlled geometric feature details to make it easy to manufacture, we developed an anisotropic MFSE model considering anisotropic spatial correlation to affect the topological complexity of the structure without requiring filters or other size-control methods. In addition, to solve the low computational efficiency arising from the existence of bound constraints in the original MFSE method in [29], a generalized sigmoid function is introduced to discard these constraints in the optimization model.

In this study, an anisotropic MFSE method is proposed to solve the optical metalens problem with the finite element method and sensitivity analysis. The anisotropic correlation lengths can be used to easily obtain simple topological features and one-dimensional solutions to avoid difficult-to-manufacture designs. Numerical examples are provided to demonstrate the validity and usability of the proposed optimization algorithm for the design of metalenses.

2. Topology representation by material field with anisotropic correlation

In continuum topology optimization, the distribution of materials is continuous and can be expressed by an uncertain material field $\varphi ({\boldsymbol x})$ with spatial correlation. To represent the structural topology, the material field is divided by $\varphi ({\boldsymbol x}) = 0$, the upper is taken as solid region ${\Omega _{\textrm{solid}}}$, and the lower is taken as void region ${\Omega _{\textrm{void}}}$, that is,

In MFSE, a correlation function is introduced to quantify the spatial correlation between any two points ${{\boldsymbol x}_i}$ and ${{\boldsymbol x}_j}$ in the design domain [29], namely,

The design domain of metalenses in this paper is a rectangular region with a relatively large aspect ratio; using the isotropic correlation length will produce many details that are difficult to manufacture in the optimized structure, so an anisotropic spatial correlation for the material field is naturally introduced;

Assume that ${N_P}$ observation points are uniformly distributed in the design domain ${\Omega _{\textrm{des}}}$. The correlation matrix ${\mathbf C}$ between any points can be constructed according to Eq.(3) as follows

The eigenvalue problem of correlation matrix ${\mathbf C}$ becomes

Through the correlation matrix ${\mathbf C}$, the material field can be expressed by a linear series expansion [31] as

Eigenvectors with very small eigenvalues contribute negligibly to the material field; thus, it is reasonable to truncate Eq.(6) to improve the computational efficiency under truncation error $\varepsilon > 0$, as follows:

Through the truncation of the material-field series, the form of $\varphi ({\boldsymbol x})$ can be reformulated using the reduced series expansion of the material field and is written as

 

Fig. 1. The first five modes of anisotropic material field.

Download Full Size | PDF

It is apparent from the first row of Fig. 1 that if the same correlation lengths in x and y directions are applied to the square domain, a circular attenuation mode is obtained in the first mode, and if the anisotropic correlation lengths are applied, the first mode will become an ellipse while the direction of the major axis of the ellipse is the direction with the larger correlation length. It can be inferred from the figure that the material field changes more slowly in a certain direction with a larger correlation length. This inspired us to use the anisotropic correlation lengths to implicitly influence the topological feature in a certain direction to facilitate manufacturing.

3. Topological design of optical metalenses

In this section, we generalize the anisotropic MFSE model based on the finite element analysis and apply it to the 2D metalens design optimization. A new projection function is introduced to avoid the bound constraints of the original MFSE method and thus improve the optimization efficiency. The topology optimization model based on the anisotropic MFSE method for metalenses design is introduced and its sensitivity is derived below.

3.1 Finite element analysis

The propagation of electromagnetic waves in the metalens is governed by Maxwell's equations. The design of the 2D case is mainly considered in this paper, while the cases from 2D to 3D are easy to derive. For the 2D case, there are two polarizations of the electromagnetic field, namely, the transverse electric (TE) and transverse magnetic (TM) modes. In the TE mode, the electric field is confined to the wave-propagation plane, and the magnetic field only has a component perpendicular to the plane. In contrast, in the TM mode, the magnetic field is confined in the wave-propagation plane, and the electric field only has a component perpendicular to the plane.

Without loss of generality, this paper mainly considers the optimization of the optical metalens under TM mode. When there is no point light source in the design domain, Maxwell's equations can be decoupled to obtain the 2D Helmholtz partial differential equations:

Commonly used numerical methods for Maxwell's equations include the plane wave-expansion method, the finite difference method, and the finite element method. Here we use the finite element method to solve the Eq. (9). To derive the finite element form, the structure needs to be discretized into several finite elements, and the electric fields can be assembled by shape functions ${\mathbf N}(r)$.

In the above formula,

By introducing the matrix ${\mathbf K}\textrm{ = }\sum\limits_e {{{\mathbf K}^e}} + {k^2}\sum\limits_e {{\varepsilon _r}(r){{\mathbf W}^e}}$, the finite element equation of the 2D metalens in TM mode can be expressed as

3.2 Anisotropic MFSE-based topology optimization model

In this work, the main purpose of designing an optical metalens is to maximize the intensity (the size of the electric field component ${{\mathbf E}_{\mathbf z}}$) at the desired position, ${r_p}$, for a given wave. The design variables are the MFSE coefficients ${\boldsymbol \eta }$ in the design domain, so the objective function is expressed in the form

The above equation is not continuous from the material field to the topology of the structure in Eq. (1), which creates numerical difficulties in the optimization process. Therefore, it is necessary to develop a smooth projection function to make this definition continuous and differentiable. Herein, to overcome the shortcomings of the original projection function used by Luo and Bao [29], which is nondifferentiable at $\varphi (x) ={\pm} 1$, we are inspired by the sigmoid function in deep learning. The sigmoid projection function is defined as

 

Fig. 2. Comparison of different projection functions. (a) The generalized sigmoid projection function; (b) The original projection function used in [29].

Download Full Size | PDF

It can be concluded from Fig. 2 that because of the discontinuity of the original projection function at $\varphi ({\boldsymbol x}) ={\pm} 1$, the material field must be bounded within $[ - 1,1]$ in the optimization model, whereas the generalized sigmoid projection function can discard the material field’s bounded constraints without numerical difficulties. This frees us from the original bound constraints of the material field, and the form of the projection function is simple enough that it does not increase computational costs. In this paper, an incremental strategy is adopted to improve numerical stability, with $\beta$ gradually increasing from 0.5 to 20.

Herein, the interpolation form of the relative permittivity is given as a linear function by the projection material field $\tilde{\varphi }({\boldsymbol x})$. Furthermore, to punish the intermediate density, an imaginary term [32] is added.

After expansion and truncation through the material field, the metalens topology optimization model based on the anisotropic MFSE method can be written as

In this optimization model, we have discarded the constraints on the bounds of the material field, which greatly improves the optimization efficiency. Furthermore, the series expansion of the anisotropic material field is solved only once before optimization, although the solution of the eigenvalue problem in Eq. (5) may take some time.

For this problem, we use the commonly used method of moving asymptotes [33] as an optimization solver, and the projection parameter is incrementally increased to ensure the stability of the optimization process stability and the final binary design. Figure 3 shows the flow chart of the anisotropic MFSE topology optimization method used to design the metalens. First, the reduced design variables are obtained after anisotropy MFSE. Then, the sensitivities are calculated based on finite element analysis and adjoint equations for the metalens. The design variables and projection parameter are updated until the objective function becomes convergent. The convergence criterion is a relative change of less than 0.5% between two subsequent generations of objective functions.

 

Fig. 3. Flow chart of the anisotropic MFSE method.

Download Full Size | PDF

3.3 Design sensitivity analysis

For a given wave, the topology optimization of the metalens requires the sensitivities of the objective and constraints relative to the design variables. In Eq. (18), the equality constraint is solved by the finite element equation, and is thus satisfied. Therefore, it is only necessary to derive the sensitivity of the objective function to the design variables by using the chain derivation rule,

The projection function in Eq. (16) yields a simple transformation:

Then, the compound function derivation rule is used to obtain

It can be seen that the derivative of the modified sigmoid is a function with itself as the dependent variable. The simplicity of its derivative is one of the reasons why it is widely used in the deep learning.

The sensitivity of the objective function to the projected material field $\frac{{\partial \,\Phi }}{{\partial \tilde{\varphi }}}$ must be separated into its real and imaginary parts, and adjoint variables are introduced for calculation.

Collecting the above formula yields that

To eliminate the first two terms in Eq. (25), which contain the derivates $\frac{{\partial {{\mathbf E}_{z,\Re }}}}{{\partial \tilde{\varphi }}}$ and $\frac{{\partial {{\mathbf E}_{z,\Im }}}}{{\partial \tilde{\varphi }}}$, the two parenthetical expressions containing adjoint variables are set to zero.

From the above formula, the equation about ${\boldsymbol \lambda }$ is

By solving the adjoint equation and substituting Eq. (27) into Eq. (24), the expression can be simplified to

4. Numerical results and discussion

In this section, numerical examples are presented to illustrate the effectiveness of the proposed optimization method. The boundary conditions of the metalens considered herein are shown in Fig. 4. Boundaries $\Gamma $ are subject to first-order absorbing boundary conditions

 

Fig. 4. Metalens design domain and boundary conditions.

Download Full Size | PDF

The relative permittivity of the solid material is set to ${\varepsilon _{r,m}} = 3$, where the relative permittivity of air is set to ${\varepsilon _{r,air}} = 1$. The equations and boundary conditions of the above problems are discretized by the finite element method based on bilinear quadratic elements, and the finite element equation can be constructed by combining the distribution of the materials in the design domain. Once focal point ${r_p}$ is determined, the field strengths are determined to obtain the objective function $\Phi $, called the figure of merit (FOM).

In Fig. 4, $\Omega $ representing rectangular domain requires finite element analysis with ${w_\Omega } \times {h_\Omega }$. We select the narrow region ${\Omega _{\textrm{des}}}$ as the design domain of metalenses; its height is ${h_{{\Omega _{des}}}}$, its width ${w_{{\Omega _{\textrm{des}}}}}$ is the same as the model area ${w_\Omega }$, and its substrate ${\Omega _S}$ is under the design domain with ${h_s} \times {w_\Omega }$.

4.1 Single-point focus problem

First, the single-point focus design problem is solved to prove the necessity and validity of the proposed method as we show in Code 1 (Ref. [34]). The structure is discretized into 400×200 finite elements with a side length of 1 nm. A substrate with ${h_s} = 20\textrm{ nm}$ is reserved at the bottom as the base of the metalens, and the height of the design domain ${h_{{\Omega _{\textrm{des}}}}} = 15\textrm{ nm}$. The focus objective point is set as ${r_p} =$ (200 nm, 120 nm), located on the central axis of the model, and the incident wavelength is 35 nm. The material field observation points are in coincidence with the centers of finite elements in the design domain, although in essence they are separated from each other.

If the isotropic correlation length in the original MFSE method is adopted, the commonly used $\alpha = 0.3$ is chosen in ${l_c} = \alpha \cdot \min ({{w_{{\Omega _{\textrm{des}}}}},{h_{{\Omega _{des}}}}} )$, which means that ${l_{cx}} = {l_{cy}} = 0.3 \times \min ({400\textrm{nm},15\textrm{nm}} )$= 4.5 nm. The corresponding optimized metalens topology is shown in Fig. 5. It can be concluded that there are isolated features in either direction due to the inadequate capabilities of the original MFSE method in handling the narrow design domain. In our numerical practice, it is still difficult to eliminate isolated topological features through the careful adjustment of $\alpha$ for this special problem; thus, the introduction of an anisotropic material field is necessary.

 

Fig. 5. Isolated topological features in the isotropic material field optimization result.

Download Full Size | PDF

When use the proposed anisotropic MFSE, the correlation length in the x-axis direction is ${l_{cx}} = 12$ nm, whereas that of y-axis direction is ${l_{cy}} = 5$ nm. The number of the MFSE design variables is 390 when the truncation error is set as ${10^{ - 4}}$, which is much smaller than the 6000-dimensional design variables of the density-based method. Figure 6(a) shows the corresponding binarization metalens configuration at this time, where black represents solid material and white represents air. Figure 6(b) indicates the final electric field distribution with the objective value of $\textrm{FOM} = 13.4740\textrm{ [}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$, which demonstrates the ability of this method to solve the optical inverse design problem. The iteration history of the objective function is given in Fig. 6(c), where it can be seen that the iteration process is robust and some slight oscillations are caused by the increase in the projection parameters.

 

Fig. 6. Single-point focus problem optimization result with ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 5 nm. (a) Metalens configuration, (b) field ${|{{{\mathbf E}_z}} |^2}$ distribution, (c) iteration history.

Download Full Size | PDF

 

Fig. 7. Single-point focus problem optimization result with ${l_{cx}} =$ 20 nm and ${l_{cy}} =$ 8 nm. (a) Metalens configuration, (b) field ${|{{{\mathbf E}_z}} |^2}$ distribution, (c) iteration history.

Download Full Size | PDF

There are some isolated features in Fig. 6(a) that are difficult to manufacture; thus, the design should be modified to meet manufacturability requirements. In the anisotropic MFSE method, the size of the correlation lengths can implicitly control the topological features. The correlation lengths are increased to ${l_{cx}} = 20$ nm and ${l_{cy}} = 8$ nm. The number of the design variables obtained in this case is 173 and the optimization process is rerun. The single-point focus problem converges to $\textrm{FOM} = 13.0060\textrm{ [}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$ after 76 iterations and yield a simpler and easier-to-manufacture structure shown in Fig. 7(a). There are no isolated features in the structure, which demonstrates the ability of anisotropic correlation length to affect the characteristics of structural topologies. To further investigate the influence of the initial design on the optimization results, two different cases of initial material-field functions are considered, namely, Case 1: $\tilde{\varphi }(x) = 0.6$; and Case 2: $\tilde{\varphi }(x) = 0.4$. The optimized topologies about the two different initial designs are shown in Figs. 8(a) and 8(b) with the final objective values of $\textrm{FOM} = 13.0993\textrm{ [}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$ and $\textrm{FOM} = 13.0050\textrm{ [}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$, respectively. It is seen that the optimized metalenses in Fig. 8 exhibit very similar structural configurations and objective performances to Fig. 7.

 

Fig. 8. Single-point focus problem optimization result with different initial designs. (a) Optimized configuration and field ${|{{{\mathbf E}_z}} |^2}$ distribution by the initial design $\tilde{\varphi }(x) = 0.6$, (b) optimized configuration and field ${|{{{\mathbf E}_z}} |^2}$ distribution by the initial design $\tilde{\varphi }(x) = 0.4$.

Download Full Size | PDF

Moreover, to verify that the material field’s bound constraints were removed via modification of the projection function, the programs were run separately on the same computer with ${l_{cx}} = 12$ nm and ${l_{cy}} = 5$ nm. The total calculation cost is 474 s with an objective value of $\textrm{FOM} = 13.5058\textrm{ [}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$ using the original projection function with the material field’s bound constraints in [29]. In contrast, the calculation cost is 291 s using the projection function presented in this paper, which indicates higher efficiency due to the removal of the material field’s bound constraints.

 

Fig. 9. Single-point focus problem optimization result that handle the visible incident wavelength of 400 nm. (a) Metalens configuration, (b) field ${|{{{\mathbf E}_z}} |^2}$ distribution, (c) iteration history.

Download Full Size | PDF

Without loss of generality, the visible light in the wavelength range of 390-780 nm can also be selected as the incident wave, but it may require a larger structural design domain. To demonstrate the ability of the proposed method in handling the visible incident wavelengths, the incident wavelength is set as 400 nm, and the entire structure is discretized into 1600×400 finite elements with the element size of 2 nm. The target point is located in ${r_p} =$ (1600 nm, 600 nm). The height of design domain and substrate are both modified to 40 nm. After rerunning the optimization process with ${l_{cx}} = 40$ nm and ${l_{cy}} = 20$ nm, the optimized result and corresponding field strength are shown in Fig. 9. The focusing effect can be seen in Fig. 9(b) without isolated features and the objective value converges to 1.9354 $\textrm{[}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}]$ as shown in Fig. 9(c). The correctness of the proposed method under the visible incident wavelengths is proved by this case.

4.2 One-dimensional design problem

In nanoscale manufacturing, the choice of manufacturing technology limits the geometric complexity of the design. For example, optical projection lithography and electron beam lithography limit the change in geometric figures to 2D patterns and then extrude them out of the plane. In the design of topological optimization of 2D optical metalenses [24], the selection of certain manufacturing technologies only enables unidirectional design, requiring additional operations in the density-based method. However, this is simple under the anisotropic MFSE method, which only needs to make the correlation length of the non-design direction much larger than the scale of the design domain.

As proof, the design domain with the width of 400 nm and the height of 15 nm is considered. We simply set the anisotropic correlation lengths as ${l_{cx}} =$ 12 nm in x-axis direction and ${l_{cy}} =$ 1500 nm in y-axis direction. Then, the first five modes of the material field series are as shown in Fig. 10(a), while the design variables are set to ${\boldsymbol \eta } = {\{{1,1,1,1,1,0,0, \cdots ,0} \}^\textrm{T}}$. The corresponding material field is plotted in Fig. 10(b) using Eq. (8). It can be seen from Fig. 10 (a) that the modes all present a clear vertical distribution, so the final material field with ${\boldsymbol \eta } = {\{{1,1,1,1,1,0,0, \cdots ,0} \}^\textrm{T}}$ will naturally show the unidirectional characteristics demonstrated in Fig. 10 (b). Thus, a simple method to solve the one-dimensional design problem is presented.

 

Fig. 10. (a) First five modes of material field with correlation lengths ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 1500 nm, (b) material field formed by the superposition of the first five material modes.

Download Full Size | PDF

To achieve the one-dimensional design, we only need to set the correlation length in y-axis direction to a large value, e.g., ${l_{cy}} =$ 1500 nm. The other parameters are consistent with the previous example in Section 4.1. The final binarized designs with ${l_{cx}} =$ 12 nm and 30 nm are shown in Fig. 11 and Fig. 12, and the objective values are 11.1129$\textrm{[}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$ and 10.2620$\textrm{[}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$, respectively. It can be seen from the optimized topologies that the design features only change in the x direction. The differences between these topologies can be discovered by comparing the optimized results. The central part in Fig. 12, with a longer correlation length, is wider than that in Fig. 11, and small features on both sides of the topology in Fig. 11 are removed in Fig. 12. Compared with the corresponding ${|{{{\mathbf E}_z}} |^2}$ field presented in Fig. 11, that in Fig. 12 confirms the focusing effect at the targeted point within the one-dimensional design problem.

 

Fig. 11. One-dimensional design problem optimization result with ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 1500 nm. (a) Metalens configuration, (b) field ${|{{{\mathbf E}_z}} |^2}$ distribution, (c) iteration history.

Download Full Size | PDF

 

Fig. 12. One-dimensional design problem optimization result with ${l_{cx}} =$ 30 nm and ${l_{cy}} =$ 1500 nm. (a) Metalens configuration, (b) field ${|{{{\mathbf E}_z}} |^2}$ distribution, (c) iteration history.

Download Full Size | PDF

4.3 Multi-point focus problem

This example shows the design of a multi-point focusing metalens, as shown in Fig. 13. The boundary conditions and incident wavelength are the same as those in the previous example, but the focus point is changed from one to multiple.

 

Fig. 13. Design domain of a multi-point focusing metalens.

Download Full Size | PDF

To optimize multi-point focusing, the average field strength of the multiple points can be maximized in Eq. (15), but the optimization of the average field strength tends to converge to a solution with very good performance at some points and poor performance at others. Therefore, the objective is defined as to maximize the worst-case performance among the points to ensure that each point reaches at least a moderate level of focusing strength, namely

The objective function of Eq. (30) is given in the form of a maximal minimum function. This function is a typical non-differentiable function and can be processed by an aggregate function. Here, we select the Kreisselmeier–Steinhauser (K-S) aggregate function [35] and transform the objective function into the following expression

The sensitivity of the modified objective function to the design variables is given in the following form:

It is worth noting that in the multi-point focusing problem, there is no additional cost of finite element calculation, that is, a single iteration still only needs to solve the finite element equation once, but the accompanying equation needs to be solved multiple times. We take five points on the central axis as an example, while the x coordinates are both 200 nm, and the $y$-coordinates are evenly distributed in an 80 nm band from 60 nm to 140 nm, i.e., ${y_i} \in \{{60\textrm{nm},80\textrm{nm},100\textrm{nm},120\textrm{nm},140\textrm{nm}} \}$. The penalty parameter of K-S aggregation is increased from 10 to 30 to ensure optimization stability. When the correlation lengths are ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 5 nm, the topology and iteration history obtained are as shown in Fig. 14. The corresponding objective function is $f = \textrm{6}\textrm{.2436[}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}]$. The electric field strengths of the five focus points from top to bottom are $\textrm{FOM} = \{{\textrm{6}\textrm{.2834,6}\textrm{.3203,6}\textrm{.3045,6}\textrm{.2849,6}\textrm{.3012}} \}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}$. It can be seen that the deviation between the field strengths of the focus points is small due to the adaption of K-S aggregate function instead of average function.

 

Fig. 14. Multi-point focus problem optimization result with ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 5 nm. (a) Metalens configuration, (b) field ${|{{{\mathbf E}_z}} |^2}$ distribution, (c) iteration history.

Download Full Size | PDF

When the correlation length is small, the optimized structure presents a more complex topology compared to the single-point focus design. Therefore, the correlation lengths are adjusted to ${l_{cx}} =$ 20 nm and ${l_{cy}} =$ 15 nm, and the corresponding optimized metalens configuration and the objective function iteration history are shown in Fig. 15. The final objective function is $f = \textrm{5}\textrm{.6256[}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$, and the electric field intensity of the five points from top to bottom is $\textrm{FOM} = \{{\textrm{5}\textrm{.6502,7}\textrm{.1387,5}\textrm{.6712,5}\textrm{.6789,5}\textrm{.7164}} \}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}$, confirming the ability of this method to solve multi-point focus design problem. In addition, through comparison of the performance of the designs in Fig. 14 and Fig. 15, the latter shows a lower ${|{{{\mathbf E}_z}} |^2}$ field strength on the concern band by the reduction of the solution space.

 

Fig. 15. Multi-point focus problem optimization result with ${l_{cx}} =$ 20 nm and ${l_{cy}} =$ 15 nm. (a) Metalens configuration, (b) field ${|{{{\mathbf E}_z}} |^2}$ distribution, (c) iteration history.

Download Full Size | PDF

4.4 Multi-wavelength focus problem

A broadband design means that the metalens can focus waves with multi-wavelength at one point, which provides better robustness in practical applications. This example demonstrates the ability of the proposed method to solve the multi-wavelength focus problem. The objective function for this problem should also be changed to keep the focusing effect as consistent as possible at each wavelength. Similar to Section 4.3, the objective function is modified to an aggregate form at multiple wavelengths. Note that a single iteration needs to solve multiple finite element and adjoint equations, which can be accelerated by parallelization. We select 30∼40 nm as the desired broadband range with five simultaneously focused wavelengths, i.e., $\lambda \in \{{30\textrm{nm},32.5\textrm{nm},35\textrm{nm},37.5\textrm{nm},40\textrm{nm}} \}$. The objective function value obtained after 70 iterations is $f = \textrm{6}\textrm{.6772[}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$ when the correlation lengths are set as ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 5 nm. The corresponding electric field strengths at the five wavelengths are $\textrm{FOM} = \{{\textrm{6}\textrm{.726,6}\textrm{.7462,7}\textrm{.0836,7}\textrm{.8979,6}\textrm{.6901}} \}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}$. Similarly, the electric field strength deviations between the five wavelengths are not large, which proves the effectiveness of the broadband design. The optimized metalens configuration is shown in Fig. 16(a) and the ${|{{{\mathbf E}_z}} |^2}$-field distributions corresponding to the five targeted wavelengths are shown in Fig. 17(b-f). It can be seen that the optimized results under all wavelengths provide a focusing effect, which confirms the correctness of the optimization result.

 

Fig. 16. Multi-wavelength focus problem optimization result with ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 5 nm. (a) Metalens configuration, (b)-(f) field ${|{{{\mathbf E}_z}} |^2}$ distribution with different incident wavelengths from 30 nm to 40 nm with $\lambda \in \{{30\textrm{nm},32.5\textrm{nm},35\textrm{nm},37.5\textrm{nm},40\textrm{nm}} \}$.

Download Full Size | PDF

 

Fig. 17. Multi-wavelength focus problem optimization result with ${l_{cx}} =$ 20 nm and ${l_{cy}} =$ 8 nm. (a) Metalens configuration, (b)-(f) field ${|{{{\mathbf E}_z}} |^2}$ distribution with different incident wavelengths from 30 nm to 40 nm with $\lambda \in \{{30\textrm{nm},32.5\textrm{nm},35\textrm{nm},37.5\textrm{nm},40\textrm{nm}} \}$.

Download Full Size | PDF

Without loss of generality, to obtain simpler features, we modify the correlation lengths to ${l_{cx}} =$ 20 nm and ${l_{cy}} =$ 8 nm and rerun the optimization process. The resulting objective function we get is $f = \textrm{6}\textrm{.4712[}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$, and the corresponding target point electric field strengths are $\textrm{FOM} = \{{\textrm{6}\textrm{.5849,6}\textrm{.5695,6}\textrm{.7944,7}\textrm{.8849,6}\textrm{.5026}} \}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}$. The corresponding metalens configuration is given in Fig. 17, presenting a good layout form that without small characteristics in which and all members of the structure are connected. The schematic diagram of the iteration history corresponding to the two results is shown in Fig. 18, where each iteration history shows a stable optimization process.

 

Fig. 18. Iteration histories with different correlation lengths, (a) Iteration history with ${l_{cx}} =$ 12 nm and ${l_{cy}} =$ 5 nm, (b) iteration history with ${l_{cx}} =$ 20 nm and ${l_{cy}} =$ 8 nm.

Download Full Size | PDF

To demonstrate that the average performance of our designed metalens over the target wavelengths is better than that for a single wavelength of 35 nm, the latter design is taken as the verification object. The incident wavelengths are swept in the interval from 25 nm to 45 nm, and compared with the optimized result for a single wavelength of 35 nm, as shown in Section 4.1. The value of the FOM in Eq. (15) as a function of wavelength is plotted in Fig. 19. It can be seen that although the optimized single-wavelength design performs best at that wavelength, the broadband metalens design performs best over the full wavelength interval, e.g., with performance greatly improved from $\textrm{FOM} = \textrm{2}\textrm{.307 [}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$ to $\textrm{FOM} = \textrm{6}\textrm{.585 [}{\textrm{V}^2}\textrm{/}{\textrm{m}^2}\textrm{]}$ at $\lambda =$ 30 nm. Moreover, this method can be extended to a wider range to achieve a broadband focus effect by accepting a performance drop at the central wavelength, if needed.

 

Fig. 19. Wavelength sweep results for broadband design and single-wavelength focus design.

Download Full Size | PDF

5. Conclusions

In this paper, the topology optimization problem of 2D optical metalenses with a narrow design domain is systematically addressed, and a new anisotropic MFSE optimization algorithm is proposed. Through finite element calculation and sensitivity analysis, the optimization algorithm focuses the electric field to a given point, and the final metalenses present no isolated features. The numerical results show that the proposed algorithm is effective and applicable to any wavelength and location. The introduction of the generalized sigmoid projection further improves the computational efficiency by discarding material field’s bound constraints in the original optimization model [29].

The main advantages of the proposed method compared to other topology optimization methods can be highlighted as follows. (1) The topological complexity of the structure can be changed easily by varying the correlation lengths of the anisotropic material field to ensure manufacturability. (2) The checkerboard patterns and mesh dependency that existed in the conventional density-based method are inherently avoided by the spatial dependency definition in the material field. (3) The anisotropic MFSE method greatly reduces the dimensionality of topological design space.

The proposed method can be easily modified to solve plasmonics inverse design problems such as metallic reflector or others by modifying the objective formulation and design sensitivity. Moreover, due to the significant reduction in the number of design variables, this method can be combined with gradient-free algorithms to investigate complex optical device design problems that are difficult to obtain design sensitivity information. In addition, an urgent challenge of nanophotonic device topological design within a narrow design domain is the structural topological integrity and its manufacturability. Numerical results show that isolated and very thin features can be avoided in optimized metalenses by varying the correlation lengths in the anisotropic material field. This will give us great facilitate to manufacture novel photonic devices by using some manufacturing technologies such as 3D laser nanolithography.