1. Introduction

In the imaging industry, there has always been an intense demand for high-performance achromatic lenses [1–17] in various applications with aberration-free imaging capability and diffraction-limited resolution. Together with the interest on lenses miniaturization from the embedded market, works focusing on achromatic ultra-thin lenses have been largely conducted [2,18–24]. Researchers also target at other parametric requirements that show great practical potential, such as large size [25,26], wide operating band [6–8,11,27], and high numerical aperture (NA) [14]. To fulfill those market hotspots, several seminal works [5,12,14,20–22,28,29] demonstrate broadband focusing, as summarized in Fig. 1(a). Some designs [12,20–22,28,29] possess large diameter $D > 500\lambda $, the extremely large achromatic lenses is also achieved with $D\sim {10^4}\lambda $ [12]. Among them, however, the achieved efficiency is not always mentioned [12,20] or not as high as expected. For high NA designs ($NA > 0.9$) with simultaneous broadband operation at 450-700 nm with an average efficiency of 23% [14], it is claimed that due to the computational capacity, the scope is confined in 2-D design with the size at tens of λ. Besides, some achromatic lenses [12,20] are optimized at three discrete wavelength points in the band instead of in continuous broadband.

 

Fig. 1. Illustration of the overview of achromatic lens design. (a) The achieved functionality of state-of-art works for achromatic flat lens, the number in the dark cycle is referred to the cited work. Bandwidth BW is defined as $({{f_u} - {f_l}} )/{f_c}$, and ${f_u},{f_l}$ and ${f_c}$ are upper, lower, and central frequency in the operating band respectively. (b) The limitation and trade-off of ultra-thin lens facing by achromatic ultra-thin lenses: metalenses (AML) [31,32], considering the group delay $\partial \psi /\partial \omega $, and ADL with both forward design (multi-order diffractive lens) [32] and inverse design [28]. $\psi $: phase produced by metaatoms, $F$: focus length, and $\mathrm{\Delta }{F_m}$: focus length different between adjacnet $m.\; {\lambda _m}:$ wavelength that MOD operates at the first order of diffraction [32], $P$: computational power, ${N_r},\; $ number of unit-zone, proportional to diameter. (c) Block diagram for our algorithm RPSA that updates the $c(\lambda )$— an array contains relative phase offset for each wavelength, during optimization. (d) The optimized achromatic lens system, with the inset plot being the normalized intensity on a line crossing the center of focal plane, and the orange arrows indicate its full width at half maximum (FWHM).

Download Full Size | PDF

It appears that there still exists vacancy for large-scale and broadband ultra-thin lenses with the acceptable efficiency defined based on the total transmitted power. This phenomenon can be analyzed separately for two categories of the achromatic ultra-thin lenses, achromatic metalenses [3,4,6,9,10,14,16,23,24,30] and achromatic diffractive lenses (ADL) [12,21,22,28,29] in Fig. 1(b). For the achromatic metalenses that explore the dispersion response from the meta-atoms with various materials and geometries, the group delay required is directly linked to the diameter D of the lenses as shown below the ‘AML’ from Fig. 1(b) [31,32]. Therefore, the achievable group delay for the meta-atom library poses direct restrictions on the operating band and diameter [6,31,32], limiting it from large-scale applications such as cameras. On the other hand, achromatic diffractive lenses based on diffractive optic elements (DOE) [18] support low-cost large-scale manufacture using lithography may provide possibility, where forward-design lenses based on multi-order diffractive lens forbids simultaneous large diameter and broadband [32] by concerning the inherent connection between the diffraction order m, multi-order focal length ${F_m}$, diffraction efficency ${\eta _d}$, and the operating wavelength $\lambda $, which is illustrated in ‘MOD’ under the ‘ADL’ in Fig. 1(b). And there is no existing analytical solution for ADL at customized wavelengths or continuously in the band. Therefore, optimization [17,33] needs to be implemented to discover the potential lens structures; but to our knowledge, the challenge for designing broadband and large-scale ADL remains unresolved. With such vacancy, we, therefore, wonder whether it is possible to achieve a high-performance ADL meeting all the industrial requirements through inverse design.

Here, we implement the inverse-design algorithm, Relative-Phase Simulated Annealing (RPSA), generalizing the process for diffraction-limited, polarization-insensitive ADL design with arbitrary sizes, operating bands, materials, and other parameters (Fig. 1(c), (d)). Compared to optimizations performed on the structure height distribution $h(r )$ [14,27], where r is the radial position, we switch the optimizing object to the relative phase offset $c(\lambda )$ for different wavelength $\lambda $. We test our RPSA on ultra-thin lenses operating in $[{8\mu m,\textrm{ }12\mu m} ]$ with diameter $D = 1cm$, exhibiting focusing spots with negligible sidelobes at all nine wavelengths optimized. However, the average focusing efficiency is numerically $9.6\%$, from which we observe that counter-intuitively, the visualized appearance of the intensity profile and numerical value of focusing efficiency are not directly correlated, especially for large-scale design. When approaching achromatic broadband, focusing efficiency dramatically drops with decreasing interval between adjacent wavelengths, and with varied design parameters, we witness the optimized focusing efficiency dropping exponentially as the diameter or bandwidth increases, where we further deduce an empirical formula to describe the pattern. Moreover, our design provides ultra-fast optimized solution especially for large-scale design, less than 10 minutes for $D$ = 10000$\lambda $ achromatic lenses at three wavelength points. We hope that our results can be used as a fast solution for further research in achromatic large-scale imaging system design, and we conclude that based on our assumptions for the algorithm, a large and broadband ultra-thin lens with high efficiency is almost impossible.

2. Design algorithm

Recently, many works for designing achromatic ultra-thin lenses based on $h(r )$ have been conducted, such as inverse design with adoint method [14], and direct binary search based on diffractive optic elements [21,22,25,27,28]. However, the large-scale optimization often encounters challenges from computational power. Focusing on spherical lenses [27] instead of cylindrical ultra-thin lenses [14] (see Appendix 6.1), we found the most time-consuming part is the calculation of the intensity profile on the focal plane, which requires one 2d Fourier transform followed by one 2d inverse Fourier transform. We tested this process for a $D\sim 1.5k{\lambda _0}$ lens in MATLAB, leading to 20000 seconds per 1000 iterations (1000 combinations of $h(r )$). Moreover, in the solution space, $h(r )$ can take ${N_r}^{N_{lvl}}$ combinations, where ${N_r}$ is the length of array r and $h(r )$, and ${N_{lvl}}$ is the number of discrete values that $h(r )$ can take. With tedious calculation and huge number of combination possibilities, large computational resource assigns hard constraints on ${N_r},$ namely, the diameter D for global optimization of large-scale ADL. To break this limitation, we instead focus on the phase profile right after the ultra-thin lenses, wherein we add a gauge term $c(\lambda )$ and engineer the figure of merit to further simplify the problem as stated in following paragraphs.

For wavelength $\lambda $ in operating range $[{{\lambda_{min}},\textrm{ }{\lambda_{max}}} ]$ and the number of wavelengths to be ${N_\lambda }$, to design a structure that provides a broadband focusing functionality, the important assumption we made is that the optimized phase profile should approximate the ideal spherical wavefront for all $\lambda $[8]:

On the other hand, the phase profile generated by the DOE structure with plane wave incidence is obtained by the phase propagation in the structure with refractive index $n(\lambda )$. For wave transmitted through the structure $h(r )$, the total optical path length at the exit plane right after the highest structure is the summation of the path in structure plus the path for exiting from structure, propagating in background media and reaching the exit plane. The phase profile is the product of optical path length and the wavevector in free space $k = 2\pi /\lambda $:

There are previous optimizations whose optimization targets is maximizing $|\cdot |^2$ [14] and $|{IFT({FT({\cdot} )} )} |^2\; $ [21,22,27,33], where $({\cdot} )$ serves as an operation of E field on exit plane, and $FT$, $IFT$ are Fourier transform and inverse Fourier transfrom when calculating the intensity on focal plane. Compared to them, our Eq. (3) is valid as all the methods reach the optimum under $\psi \; = {\psi _{ideal}}$. The reason for choosing optimizing $Re({\cdot} )$ is for the convenience of coding. In addition, by replacing the non-linear term in optimization target with $Re({\cdot} )$, the problem is greatly simplified because the structure at each r can be optimized individually as elaborated in discussion section.

Here we implement a Simulated Annealing (SA) in MATLAB for global optimization of $c(\lambda )$. Compared with the local optimization algorithms such as the gradient-based inverse design with adjoint method [14], and direct binary search [33], the randomness included in SA guarantees a nearer globally optimized result. Our SA optimization operates on $c(\lambda )$ which is a relative phase term in the Eq. (1) and thus called Relative-Phase Simulated Annealing (RPSA). Briefly shown in Fig. 1(c), all the design parameters (in Eq. (1)–(3)) are provided to the algorithm and the relative phase $c(\lambda )$ is randomly initated. When iterating, for each $c(\lambda ),$ the $h(r )$ with individual ${r_i}$ can be optimized independently for an overall maximum of Eq. (3). At each ${r_i}$, the optimized $h({{r_i}} )$ from ${N_{lvl}}$ values can be directly calculated using the method of exhaustion. By combining ${r_i}$, the optimized $h(r )$ for this $c(\lambda )$ is obtained, where the $FOM$ can be calculated based on Eq. (4). From the SA perspective, at iteration m with updating the relative phase ${c_m}(\lambda )$, the optimized ${h_m}(r )$ is also generated with $FO{M_m}$. The SA accepts the ${c_m}(\lambda )$ based on both the $\mathrm{\Delta }FOM = FO{M_m} - FO{M_{m - 1}}$ and the current iteration number m (related to the term of temperature in simulated annealing algorithm) which follows the conventional SA flow [34]. The RPSA iteration ends when converging or reaching the maximum iteration number. The RPSA will generate the last $c(\lambda )$ with its associated $h(r )$, which will be the final optimized structrure for our ADL.

For the running time, in contrast to 20000 seconds, for the same $D\sim 1.5k{\lambda _0}$, the RPSA takes only ${\sim} 5$ seconds for 1000 iterations. Together with a smaller solution space, $c(\lambda )$ as variables compared to $h(r )$, which also decreases the number of iterations required, RPSA bypasses the challenge of computational requirement for large-scale optimization. Moreover, our approach can be applied for smaller lenses, as long as the assumptions for DOE are valid, providing a new platform for ADL with custom configurations without size constraints. Generally, the design requirements such as geometry, material, and frequency, and specifications for algorithms such as converging speed can be customized. One thing should be noted is that our solution space based on $c(\lambda )$ is a subset of the one based on $h(r )$, which is confined by Eq. (3), and therefore, we admit the possibility for an optimization tool such as neural network to find a better solution outside our solution space which do not obey Eq. (3). However, due to the huge computation demand for an ${N_r}$-dimension solution space, the directly searching for $h(r )$ tends to converge at a very localized solution, and to our knowledge, there is still no work present an achromatic broadband large-scale ultra-thin lens with high efficiency.

The criterion we use for evaluation is the “focusing efficiency”, defined as the intensity ratio of the total power within the main lobe at the focal plane to the total transmitted power through the structure $h(r )$. We want to note that as a diffraction-limited lens, for numerator, the diameter for the main lobe is approximately three times of full width at half maximum (FWHM), indicating our focusing efficiency and those defined based on $3 \times FWHM\; $[15,22] should share similar values numerically. Moreover, unlike the works whose focusing efficiency has a reference denominator power as the total incident power [25,26] onto the sensor from the measurement (sometimes the diameter is $({10 + } )\times FWHM$ only) [35], leading to a high numerical value for the “focusing efficiency”, our definition based on total transmission power (integrated on $({ > 10^3} )\times $ FWHM for large scale lens) possess a more physical meaning to evaluate the overall performance for the system.

3. Results

3.1 Large-scale broadband ADL optimized on nine wavelength points

To study the optimization performance of RPSA, first, we design an achromatic large-area lens operating in LWIR of $8\mu m - 12\mu m,$ with thickness of $10\mu m$, diameter $D = 1000{\lambda _0}$, where ${\lambda _0} = 10\mu m$ is the central wavelength. As shown in Fig. 2(a), we set the unit size of the structure $h(r )$ in radial direction as $10\mu m$, resulting in $500$ zones. The $NA$ is $0.3714$, and thus focal length $F = 12.5mm$. The material is silicon with $n(\lambda )= 3.42$ over the operating band. The optimization is performed on $9$ uniformed distributed wavelengths ${\lambda _i},\; i = 1,2, \ldots ,9$.

 

Fig. 2. Performance of the achromatic lens design with diameter of $1000{\lambda _0}$ and optimized at nine wavelength points. (a) The layout of the optimized structure. (b) The optimized relative phase $c(\lambda )$ for nine wavelength points. (c) The optimized intensity profiles at focal plane. (d) The simulated intensity profiles using Lumerical FDTD. (e) The optimized and simulated FWHMs and focusing efficiency.

Download Full Size | PDF

The optimized structure from RPSA is shown in Fig. 2(a), with the optimized $c(\lambda )$ in Fig. 2(b). Here we let the $c(\lambda )$ be selected from $30$ numbers uniformly distributed in $[{ - \pi ,\pi } ]$. $c(\lambda )$ and $\lambda $ present a monotonic increasing relation. To analyze the intensity on the focal plane, its electric field distribution $E({x,y,z = F} )$ is calculated through the Rayleigh-Sommerfeld theory [36], while the field on exit plane $E({x,y,z = 0})$ is attained from Eq. (2).

The intensity at $z = F$ is presented in Fig. 2(c) and 2(d). At all optimized ${\lambda _i}$, the structure provides a diffraction-limited focusing spot with negligible side lobes. The FWHMs plotted as blue circle dashed line in Fig. 2(e) showing near-diffraction-limited main lobe width. Despite the good appearances, the average focusing efficiency, our criteria for the evaluation, is relatively low $:9.6\%$, which is the mean value of the optimized efficiency in Fig. 2(e). This low efficiency is intuitively expected as the design target is challenging with both large area and broadband, and can be explained as although the area outside the main lobe is dark (low intensity) in Fig. 2(c), the lens size is ${10^6}\lambda _0^2$, whereas for the main lobe area is only around $10\lambda _0^2$. The denominator part of focusing efficiency calculated from integrating the intensity (low value) over the entire area ${\sim} {10^6}\lambda _0^2$ may result in high value, and thus causes a low average focusing efficiency. Moreover, if we comprise our design target by narrowing down the band or reducing device size, a much higher efficiency should be expected.

Before the discussion of efficiency dependence of the design parameters, we first validate this design through simulation in Lumerical FDTD. As shown in Fig. 2(e) with the solid line, the FDTD simulated focusing efficiency agrees with the optimized focusing efficiency, indicated with the dashed line, as is the FWHM. Comparing Fig. 2(c) and Fig. 2(d), we can conclude that our assumptions for the algorithm are acceptable. In addition, from the simulation, the structure shows a transmission efficiency of 70%, leading to an overall efficiency of $70\%\ast 9.6\%= 6.72\%$.

3.2 Varying parameters and focusing efficiency trends

To study the capability of our algorithm, we vary the design parameters to test its performance for discrete-wavelength-based ADL. Keeping the other parameters unchanged, we modify the number of units ${N_r}$, thus the diameter D as well as the bandwidth $\mathrm{\Delta }\lambda $. Remarkably, when dealing with non-ideal focusing profiles, it is hard to define the ‘main lobe’. Therefore, in the following sections, the focusing efficiency is re-defined as the intensity within the diameter of 3 times of the ideal FWHM as discussed before, with $FWH{M_{ideal}} = \lambda /2NA$, divided by the total transmitted power, eliminating the disturbance of main lobe width and the non-ideal shape such as the concave at the peak.

Intuitively, it’s easier to design a focusing device with a small size or narrow bandwidth, resulting in higher focusing efficiency. Together with our inverse-design results with different parameters, the optimized performance can be roughly divided into three sections according to the difficulty, and we just name them as ‘easy, ‘medium’, and ‘challenging’ cases. For the easy case, the diameter D and/or bandwidth $\mathrm{\Delta }\lambda $ are very small, e.g., $({40\; {\lambda_0},\; 4\mu m} )$ (shown in Fig. 3(a)) and $({1000\; {\lambda_0},\; 0.02\mu m} )$. In these cases, the algorithm finds the required structure easily, achieving near diffraction-limit focusing performance with high focusing efficiency $({ > 50\%} )$. Those structures usually present continuous broadband focusing functionality; namely if we interpolate 9 ${\lambda _i}$ by adding m point between each adjacent wavelength pair and test at $({9 + ({9 - 1} )\times m} )$ points, all data points exhibit similar performance. Bottom panel of Fig. 3(a) demonstrates the $m = 1$ case, with the interpolated $\lambda $ marked with filled markers in different color, showing near diffraction-limit focusing performance. The reason for the optimized FWHM being slightly higher than the ideal value $({\lambda /2NA} )$ is that, as shown in Fig. 2(a), the height distribution is discrete rather than continuous in both radial direction (discrete zones) and $z - \; $ direction (discrete height level).

 

Fig. 3. PSF, focusing efficiency and FWHM of the optimized lenses with different parameters, tested at interpolation points. The optimized wavelength points ${\lambda _i}$ are marked as the title for sub-plots, while no title for plots at interpolated points. For scenarios when we cannot define a main lobe and therefore FWHM, the FWHM value is not marked in plots in bottom panel. (a) for the easy case with either small diameter or bandwidth, e.g. $D = 40{\lambda _0},\; \; \mathrm{\Delta }\lambda = 4\mu m$. (b) for the challenging case with large diameter and bandwidth, e.g. $D = 1000{\lambda _0},\; \; \mathrm{\Delta }\lambda = 4\mu m$, and (c) for the medium case with medium diameter and/or bandwidth, e.g. $D = 1000{\lambda _0},\; \; \mathrm{\Delta }\lambda = 0.8\mu m$.

Download Full Size | PDF

For the challenging case, such as the example in Fig. 3(b), the diameter and/or bandwidth have large values, e.g., $({1000\; {\lambda_0},\; 4\mu m} )$, leading to a hard achieving continuously broadband and low-efficiency performance. Meanwhile, with discrete wavelength points to optimize, the wavelength interval between ${\lambda _i}\; $ is relatively wide. Therefore, instead of a broadband system, the algorithm automatically treats the problem as a multi-functional focusing lens with different wavelengths. Consequently, for the optimized points, the focus is formed with diffraction-limit resolution, while for the interpolated points, the efficiencies are very low, with non-Gaussian the intensity profiles.

For the medium case, the diameter and/or bandwidth are larger than the easy case and smaller than the challenging case, $({1000\; {\lambda_0},\; 0.8\mu m} )$ in Fig. 3(c). As the algorithm finds it challenging to have a continuous broadband focusing and the separations between each operating wavelength are not large enough to be treated independently, RPSA makes a tradeoff between the continuously broadband (easy case) and the performance of the discrete points (challenging case). As a result, the shapes of the intensity profiles for optimized are not as good as the former two cases and the interpolated points show similar but slightly worse performance compared to the adjacent optimized points.

Keeping other parameters constant, the optimized focusing efficiency with respect to diameter (constant $\mathrm{\Delta }\lambda = 4\mu m$) and bandwidth (constant $D = 1000\; {\lambda _0}$) are shown in Fig. 4(a) and Fig. 4(b), respectively. The blue lines in both plots show an exponentially decreasing relationship between the optimized focusing and the varied parameter. We mark the focusing efficiency $\eta = 20\%$ as the horizontal yellow dashed lines. As a result, we may state that in this multi-level diffractive element approach, a large area $({ > 1000 {\lambda_0}} )$ and broadband lens with high focusing efficiency is not achievable, even considering the free shift of $c(\lambda )$. Moreover, we want to emphasize that some optimization results suffer from interpolation issues, namely, except for the optimized ${\lambda _i}$, the interpolated points present bad focusing performance. This interpolation issue becomes severer when the diameter and/or the bandwidth increases. The overall focusing efficiencies, tested at $({9 + ({9 - 1} )\times 3} )\; $ points are also plotted in Fig. 4(a) and Fig. 4(b) in the red curve and are in accordance with the previous analysis. If we term the circumstance with strong interpolation issue in Fig. 3(b) as ‘discrete broadband’, in comparison, the case in Fig. 3(a) will be ‘continuous broadband’ without interpolation issue. Moreover, the numbers of optimized points ${N_\lambda }$ also play roles in the final performance as shown in Fig. 4(c), which states that increasing the ${N_\lambda }$ leads to a continuously broadband performance (optimized efficiency approaches optimized + interpolated efficiency).

 

Fig. 4. Focusing efficiency drops exponentially with varying parameters, for both discrete (a-c) and continuous (d-e) broadband using RPSA. (a-b), With number of discrete wavelength points ${N_\lambda } = 9$, the focusing efficiency vs varying diameter D and bandwidth $\mathrm{\Delta }\lambda $ respectively. The blue lines represent the efficiency at the optimized wavelength points while the red lines also includes the interpolation points and same for (c), the focusing efficiency with varied ${N_\lambda }$. (d-e), Empirical formula fitted for several cases using exponential equations for continuous broadband. The fitted curves are shown as solid lines while the data as triangle markers. Surprisingly, the focusing efficiency optimized as continuous broadband in both plots can be merged into one fitted empirical formula: $\eta = \textrm{exp}({ - 0.073\; {D^{1/2}}\mathrm{\;\ \Delta }{\lambda^{1/3}}} )$.

Download Full Size | PDF

3.3 Broadband empirical formula for focusing efficiency of ADL

The fundamental difficulty encountered when realizing broadband system is that from the same focusing structure, different wavelengths converge at different focuses along propagation direction, leading to chromatic aberration [37]. For achromatic lenses design, exploring the extended depth of focus [13,19] with non-spherical wavefront achieves continuous broadband at the expense of focusing performance (Appendix 6.3). For our discussion on continuous broadband, we stick to the spherical wavefront in Eq. (1). As mentioned before in Fig. 4(c), increasing number of wavelength points to optimize $({{N_\lambda }} )$ offers opportunity to achieve near continuously achromatic, which is the methodology we adopt in this section.

We assign ${N_\lambda }$ with different values for different parameter setups to achieve continuous broadband. The two plots for the focusing efficiency with respect to varied diameter and bandwidth, respectively, are shown in Fig. 4(d) and 4(e). All the plots reveal similar trends with the discrete type as in Fig. 4(a–c). For Fig. 4(d), when varying the diameter, the bandwidth is fixed as $1\mu m$, and $4\mu m$, and for Fig. 4(e), when varying the bandwidth, the diameter is $200\; {\lambda _0}$, and $400\; {\lambda _0}$. In Fig. 4(d), we choose the exponential equation as the curve to fit the data points.

Surprisingly, when varying the diameter in Fig. 4(d), the two cases both result in a value of ${\sim} 1$ for the parameter a, which has the physical meaning that as the diameter reduces, the design target becomes simpler, and the achievable efficiency will finally converge to 100% as ideal. This intriguing physical explanation provides larger possibility for our fitting model to be intuitive and realistic.

To merge bandwidth, diameter, and focusing efficiency into a single equation, with Eq. (4) and $a = 1$, we fit Fig. 4(e) with

To further validate the fitting curve, for Fig. 4(d), $\mathrm{\Delta }\lambda = 1\mu m$ and $4\mu m$, surprisingly we find the ratio between b shown in the legend of the figure is $\frac{{0.11}}{{0.07}}\sim 1.6 \approx {\left( {\frac{{4\mu m}}{{1\mu m}}} \right)^{\frac{1}{3}}}$, where $(\; )^1/3\; $ is the fitted value for d, and for Fig. 4(e), $D = 200\textrm{ }{\lambda _0}$, and $400\textrm{ }{\lambda _0}$, the ratio between c shown in the legend is $\frac{{1.43}}{{1.04}}\sim 1.4 \approx {\left( {\frac{{400{\lambda_0}}}{{200{\lambda_0}}}} \right)^{\frac{1}{2}}}$, where $(\; )^1/2$ is in Eq. (4). Therefore, from Eq. (4) and Eq. (5), we may conclude with an empirical formula for focusing efficiency:

4. Discussion

For the computational challenge in large-scale lenses design, one is that time-usage explodes for larger $D\; $ when calculating the focal intensity based on the E field on exit plane, with complexity $O({{N_r}^2log({N_r^2} )} )$ for 2D Fourier transform in each epoch, and the other is that dimension of solution space increases with the increased D, both solved using RPSA. First, in Eq. (3), the summations with respect to r and $\lambda $ are linear operators. Namely, given fixed $c(\lambda )$, the optimized $h(r )$ from ${N_{lvl}}$ values can be directly calculated using the method of exhaustion at each r as we shown in Code 1 (Ref. [38]). In other words, fixed $c(\lambda )$ will generate fixed $h(r )$ without calculating the corresponding electrical field and intensity. Consequently, we bypass the intensity calculation with $O({{N_r}^2log({N_r^2} )} )$, instead, applying the method of exhaustion whose computational-usage increases linearly with respect ${N_r}(D )$ (Fig. 5 in blue dash line), which is much more time-efficient.

 

Fig. 5. Time usage of RPSA for large-scale design with respect to diameter. The orange line represents the total optimization time used and the blue line represents the average time per epoch, and both show a near-linear relationship which is controllable.

Download Full Size | PDF

Secondly, in our approach, the variables to optimize have shifted from $h(r )\; $[14,33] to the gauge term $c(\lambda )$, indicating that the dimension of solution space implied by the number of variables to optimize is also changed from ${N_r}\; $ to ${N_\lambda },$ which is independent of the lens size. For large-scale lens, commonly ${N_\lambda } \ll {N_r}$, thus our approach significantly simplifies the problem handled over to the algorithm, breaking the constraint on the lens size r, reducing the computation time and increasing the possibility for achieving global optimization. Therefore, compared with optimization on $h(r )$, not only the computation speed but also the focusing efficiency are enhanced for the same design specification. Demonstrating the large-scale optimization capability of our RPSA, in Fig. 5 with orange dash line, the ultra-large $D = 10000{\lambda _0}$ lenses with ${N_\lambda } = 3$ takes only several minutes to optimize, showing our approach offers simple and time-efficient method and platform for arbitary ADL discussions (Appendix 6.4–6.5).

5. Conclusion

To summarize, our RPSA provides a fast tool for the most diffractive achromatic lenses design with arbitrary specifications, especially helpful in large-scale design. Moreover, taking advantage of the small amount of computational source required, the relationship between focusing efficiency and bandwidth and diameter is deduced with an empirical formula, which pads the not-fully explored inverse-design ADL. The presented exponentially dropping pattern may give hints to the vacancy of large-scale broadband ultra-thin lenses in the market. We hope that the methodology of RPSA to speed up the optimization can be utilized to other optimization problems and those findings of focusing efficiency of ultra-thin lenses can provide guidance and assistance in future research and various applications for achromatic lenses industry.

6. Appendix

6.1. Cylindrical lens and spherical lens

Here we want to discuss the possibility of different kinds of lenses designed based on the optimized one-dimension array $h(r )$, where r is the radius. Our optimization target is to push the phase produced by the structure closer to the ideal spherical phase distribution with the single dimension of r regardless of the other dimensions. Therefore, the optimized $h(r )$ can be used to produce both cylindrical lenses, $h({x,y} ),$ invariant in $y$-direction, and spherical lenses that are used in previous sections $h({r,\theta } )$, invariant in $\theta $-direction.

To verify this concept, with multiple sets of design parameters, we modify the optimized results $h(r )$ with a random error term $e(r )$ with different weights and observe the focusing efficiency (within $3FWH{M_{ideal}}$) for both cylindrical lenses and spherical lenses. The efficiency in Fig. 6 shows a monotonic increasing relationship between spherical lenses efficiency and cylindrical lenses efficiency, demonstrating that a good rotational-symmetry structure for focusing usually offers good cylindrical focusing as well by modifying the structure from $\theta $-invariant $h({r,\theta } )$ to $y$-invariant $h({x,y} )$.

 

Fig. 6. The performance of a cylindrical lens and spherical lens based on same optimized structure presents monotone increasing relationship.

Download Full Size | PDF

6.2. Some details about Eq. (3)

Consider cylindrical lens case with x only, first we choose to use the following optimization target: minimizing the difference between the optimized E field $E({x,FL} )$ and the ideal diffraction limited E field ${E_{ideal}}({x,FL} )$ at the focal plane ($z = FL$):

$$\mathop \smallint \limits_{ - \infty }^\infty {|{E({x,FL} )- {E_{ideal}}({x,FL} )} |^2}dx\; $$

$$= \mathop \smallint \limits_{ - \infty }^\infty {|{E({x,FL} )} |^2} + {|{{E_{ideal}}({x,FL} )} |^2} - 2\; Re({{E^\ast }({x,FL} ){E_{ideal}}({x,FL} )} )dx\; $$

$$= \mathop \smallint \limits_{ - \infty }^\infty {|{E({x,FL} )} |^2}dx + \mathop \smallint \limits_{ - \infty }^\infty {|{{E_{ideal}}({x,FL} )} |^2}dx - \mathop \smallint \limits_{ - \infty }^\infty 2\; Re({{E^\ast }({x,FL} ){E_{ideal}}({x,FL} )} )dx$$

where the first two terms are linked to the total intensity and are with fixed value. Therefore, the design target has been simplified to maximize the $Re\left( {\mathop \smallint \limits_{ - \infty }^\infty {E^\ast }({x,FL} ){E_{ideal}}({x,FL} )dx} \right)$

Consider the overlap $\xi ({z = FL} )$ between the optimized E and the ideal target E at focal plane:

$$\xi ({z = FL} )= \mathop \smallint \limits_{ - \infty }^\infty {E^\ast }({x,FL} ){E_{ideal}}({x,FL} )dx$$

We can prove that $\xi (z )= \; \xi ({z = 0} )$ using Parseval's theorem and Rayleigh-Sommerfeld theory with the help of the impulse-response function. As a result, we only need to consider the field at $z = 0$, namely, the exit plane (Eq. (2)). Here the optimization target has been reduced to

$$Re\left( {\mathop \smallint \limits_{ - \infty }^\infty {E^\ast }({x,0} ){E_{ideal}}({x,0} )dx} \right)$$

For the algorithm which is discrete in x, $\mathop \smallint \limits_{ - \infty }^\infty {E^\ast }({x,0} ){E_{ideal}}({x,0} )dx$ can be written as$\mathop \sum \limits_x \textrm{exp}({ - i\textrm{ }{\psi_{ideal}}(x )} )\times \textrm{exp}({i\textrm{ }\psi (x )} )$, where $\psi $ means the phase at exit plane.

Consider multi-wavelength also, maximize $\xi ({z = 0} )$ is maximize

$$Re\left( {\mathop \sum \limits_{x,\lambda } \textrm{exp}({ - i\textrm{ }{\psi_{ideal}}({x,\lambda } )} )\times \textrm{exp}({i\textrm{ }\psi ({x,\lambda } )} )} \right)$$

We argue that by considering the Appendix 6.1, a good cylindrical performance will lead to a good spherical performance. Therefore, the Eq. (3) from the main text is derived with the x changed to r.

6.3. Extended depth of focus (EDOF) and continuously broadband

Another interesting functionality for focusing-related devices is extended depth of focus (EDOF). We observe the pseudo-dyadic relation between the EDOF and continuously broadband system. In other words, the structure EDOF always demonstrates a possessing continuously broadband focusing functionality at single wavelength point in the middle of the operating band as presented in Fig. 7(a). Moreover, for a focusing lens at several discrete wavelengths as shown in Fig. 7(b), the $x - z$ plot displays its axial performance, showing a dash line near the focal point as this discrete focus and the discrete wavelength points also share the pseudo-dyadic relation. Moreover, in Fig. 8(a), compared to the PSF that optimized at nine wavelength points with RPSA in Fig. 3(b), the PSF for the design based on EDOF [13] with same design parameters in Fig. 8(b) presents a continuously broadband performance. Consequently, this observation offers another approach to achieve a continuously broadband focusing designs by designing the EDOF structures with non-spherical phase profiles [13,19]. In those approaches, their main lobes at the intensity profiles keep dominant compared to the sidelobes, sacrificing the full width at half maximum (FWHM). Moreover, from Fig. 8, the EDOF approach designs also face the challenge of lower focusing efficiency compared to RPSA approach for the same setup at large diameter and bandwidth.

 

Fig. 7. The axial intensity profiles towards z direction shows (a) for EDOF, the axial intensity profile presents as a continuous focal (optimized at 9 points), (b) discrete focal points for discrete wavelength optimization using RPSA.

Download Full Size | PDF

 

Fig. 8. Normalized intensity on focal plane, focusing efficiency and FWHM of the optimized lenses with different parameters, tested at interpolation points. For scenarios when we cannot define a main lobe and therefore FWHM, the FWHM value is not marked in plots in bottom panel. (a) same as Fig. 3(b) in main text, optimized points show focal performance and the interpolated points show poor efficiency. (b) for the EDOF, all wavelength presents focal performance, but shows poor FWHM and low efficiency.

Download Full Size | PDF

6.4. Comparison between different optimization methods

For designing an achromatic lens with the help of the algorithm, there are works before using $|\cdot |^2$ at exit plane [14] and intensity of $IFT({FT({\cdot} )} )$ representing the propagation towards the focal plane [33] as the optimization target to ignore the phase offset at the focus, where $({\cdot} )$ serves as an operating of E field on exit plane right after the structure. However, for each iteration, the Fourier Transforms takes a long time and especially for spherical lenses that considering both $({r,\; \theta } )$. When the diameter $D$ becomes larger, this timing issue becomes severer, and even for $|\cdot |^2$ method that ignores the propagation, the number of possible solutions in solution space that depends on $D$ becomes huge $({N_r^{{N_{lvl}}}} )$, resulting in the algorithm difficult to converge and tending to find a local optimized solution even adopting a global optimization method.

Instead of using $|\cdot |^2,$ we use $\mathop \sum \limits_{r,\lambda \; } Re({\cdot} )$ as the optimization target. As $\mathop \sum \limits_{r,\lambda \; } Re({\cdot} )$ are all linear operators and their order can be switched during calculation, each ${r_i},$ $({{i_{max}} = {N_r} = 1000} )$ in our example in Fig. 2 can be treated independently, and the problem is separated into i sub-problems where each subproblem has a small solution space of $1 \times {N_{lvl}}$. Those simple sub-problems can be directly solved by the method of exhaustion and thus greatly reduce the required time. In addition, the freedom of phase at the focal plane is handled by an additional factor $c(\lambda )$ as described in the method sections. Compared to $h(r )$ used in most cases and especially for large achromatic lenses design, this $c(\lambda )$ has a much smaller variable space ${N_\lambda }$ than ${N_r}$. As a result, our approach converges much faster and provides a more global result.

Using the same global optimization method (simulated annealing), we try to optimize an achromatic lens with a diameter of 1520 ${\lambda _0}$, a band of $({8\mu m,\; 12\mu m} )$, and at 9 discrete wavelengths points. The time to test $5 \times {10^4}\; $ structures is $25 \times {10^4}s$ for ${|{IFT({FT({\cdot} )} )} |^2}\; $ method, $90s$ for $|\cdot |^2$ method, and the resulting average focusing efficiency is $2.7\%$ and $3.5\%$ respectively. For our algorithm, the height distribution $h(r )\; $ is directly calculated, and the time for testing different structures is not required, achieving a $7.1\%$ focusing efficiency, around twice the value of the other two methods.

6.5. Simulation result for${\boldsymbol D} = {10^5}{{\boldsymbol \lambda }_0}$ lens

The performance of optimized lens in section 4 of main text with ${\boldsymbol D} = {10^5}{{\boldsymbol \lambda }_0}$ and ${{\boldsymbol N}_{\boldsymbol \lambda }} = 3{ }$ is shown in Fig. 9. To balance the accuracy and the huge computation power required for calculating the intensity profile of the ${\boldsymbol D} = {10^5}{{\boldsymbol \lambda }_0}$ structure, the resolution used for calculation is very low. Nevertheless, the focal performance with good FWHM (numerically value not accurate with low calculation accuracy) and efficiency is presented.

 

Fig. 9. Normalized intensity at focal plane, focusing efficiency and FWHM of the optimized lenses with $D = {10^5}{\lambda _0}$, optimized at three wavelength points.

Download Full Size | PDF

Funding

Agency for Science, Technology and Research (A2083c0060).

Acknowledgement

We acknowledge C. Z and G. H for the helpful discussion.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

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. Y. Li, X. Li, M. Pu, Z. Zhao, X. Ma, Y. Wang, and X. Luo, “Achromatic flat optical components via compensation between structure and material dispersions,” Sci. Rep. 6(1), 1–7 (2016). [CrossRef]  

2. M. Ye, V. Ray, and Y. S. Yi, “Achromatic Flat Subwavelength Grating Lens over Whole Visible Bandwidths,” IEEE Photonics Technol. Lett. 30(10), 955–958 (2018). [CrossRef]  

3. R. J. Lin, V. C. Su, S. Wang, M. K. Chen, T. L. Chung, Y. H. Chen, H. Y. Kuo, J. W. Chen, J. Chen, Y. T. Huang, J. H. Wang, C. H. Chu, P. C. Wu, T. Li, Z. Wang, S. Zhu, and D. P. Tsai, “Achromatic metalens array for full-colour light-field imaging,” Nat. Nanotechnol. 14(3), 227–231 (2019). [CrossRef]  

4. Z. Shi, M. Khorasaninejad, A. Y. Zhu, W. T. Chen, V. Sanjeev, A. Zaidi, and F. Capasso, “Achromatic metalens over 60 nm bandwidth in the visible,” in 2017 Conference on Lasers and Electro-Optics, CLEO 2017 – Proceedings (2017), pp. 1–2.

5. M. Khorasaninejad, F. Aieta, P. Kanhaiya, M. A. Kats, P. Genevet, D. Rousso, and F. Capasso, “Achromatic Metasurface Lens at Telecommunication Wavelengths,” Nano Lett. 15(8), 5358–5362 (2015). [CrossRef]  

6. S. Shrestha, A. C. Overvig, M. Lu, A. Stein, and N. Yu, “Broadband achromatic dielectric metalenses,” Light: Sci. Appl. 7(1), 85 (2018). [CrossRef]  

7. Z. B. Fan, H. Y. Qiu, H. L. Zhang, X. N. Pang, L. D. Zhou, L. Liu, H. Ren, Q. H. Wang, and J. W. Dong, “A broadband achromatic metalens array for integral imaging in the visible,” Light: Sci. Appl. 8(1), 1–10 (2019). [CrossRef]  

8. W. T. Chen, A. Y. Zhu, V. Sanjeev, M. Khorasaninejad, Z. Shi, E. Lee, and F. Capasso, “A broadband achromatic metalens for focusing and imaging in the visible,” Nat. Nanotechnol. 13(3), 220–226 (2018). [CrossRef]  

9. Q. Cheng, M. Ma, D. Yu, Z. Shen, J. Xie, J. Wang, N. Xu, H. Guo, W. Hu, S. Wang, T. Li, and S. Zhuang, “Broadband achromatic metalens in terahertz regime,” Sci. Bull. 64(20), 1525–1531 (2019). [CrossRef]  

10. H. Zhou, L. Chen, F. Shen, K. Guo, and Z. Guo, “Broadband achromatic metalens in the midinfrared range,” Phys. Rev. Appl. 11(2), 024066 (2019). [CrossRef]  

11. S. Wang, P. C. Wu, V. C. Su, Y. C. Lai, C. Hung Chu, J. W. Chen, S. H. Lu, J. Chen, B. Xu, C. H. Kuan, T. Li, S. Zhu, and D. P. Tsai, “Broadband achromatic optical metasurface devices,” Nat. Commun. 8(1), 1–9 (2017). [CrossRef]  

12. S. Wang, C. Zhou, Z. Liu, and H. Li, “Design and analysis of broadband diffractive optical element for achromatic focusing,” in Proceedings of SPIE - The International Society for Optical Engineering(2017).

13. L. Huang, J. Whitehead, S. Colburn, and A. Majumdar, “Design and analysis of extended depth of focus metalenses for achromatic computational imaging,” Photonics Res. 8(10), 1613–1623 (2020). [CrossRef]  

14. H. Chung and O. D. Miller, “High-NA achromatic metalenses by inverse Design,” Opt. Express 28(5), 6945–6965 (2020). [CrossRef]  

15. F. Balli, M. Sultan, S. K. Lami, and J. T. Hastings, “A hybrid achromatic metalens,” Nat. Commun. 11(1), 3892 (2020). [CrossRef]  

16. F. Zhao, X. Jiang, S. Li, H. Chen, G. Liang, Z. Wen, Z. Zhang, and G. Chen, “Optimization-free approach for broadband achromatic metalens of high-numerical-aperture with high-index dielectric metasurface,” J. Phys. D: Appl. Phys. 52(50), 505110 (2019). [CrossRef]  

17. J. Hu, C. H. Liu, X. Ren, L. J. Lauhon, and T. W. Odom, “Plasmonic Lattice Lenses for Multiwavelength Achromatic Focusing,” ACS Nano 10(11), 10275–10282 (2016). [CrossRef]  

18. T. Grulois, G. Druart, N. Guérineau, A. Crastes, H. Sauer, and P. Chavel, “Extra-thin infrared camera for low-cost surveillance applications,” Opt. Lett. 39(11), 3169–3172 (2014). [CrossRef]  

19. L. Zhang, J. Ding, H. Zheng, S. An, H. Lin, B. Zheng, Q. Du, G. Yin, J. Michon, Y. Zhang, Z. Fang, M. Y. Shalaginov, L. Deng, T. Gu, H. Zhang, and J. Hu, “Ultra-thin high-efficiency mid-infrared transmissive Huygens meta-optics,” Nat. Commun. 9(1), 1–9 (2018). [CrossRef]  

20. H. C. Wang, C. H. Chu, P. C. Wu, H. H. Hsiao, H. J. Wu, J. W. Chen, W. H. Lee, Y. C. Lai, Y. W. Huang, M. L. Tseng, S. W. Chang, and D. P. Tsai, “Ultrathin Planar Cavity Metasurfaces,” Small 14(17), 1703920 (2018). [CrossRef]  

21. M. Meem, S. Banerji, A. Majumder, F. G. Vasquez, B. Sensale-Rodriguez, and R. Menon, “Broadband lightweight flat lenses for long-wave infrared imaging,” Proc. Natl. Acad. Sci. U. S. A. 116(43), 21375–21378 (2019). [CrossRef]  

22. S. Banerji, M. Meem, A. Majumder, C. Dvonch, B. Sensale-Rodriguez, and R. Menon, “Single flat lens enabling imaging in the short-wave infra-red (SWIR) band,” OSA Continuum 2(10), 2968–2974 (2019). [CrossRef]  

23. H. Zuo, D. Y. Choi, X. Gai, P. Ma, L. Xu, D. N. Neshev, B. Zhang, and B. Luther-Davies, “High-Efficiency All-Dielectric Metalenses for Mid-Infrared Imaging,” Adv. Opt. Mater. 5(23), 1700585 (2017). [CrossRef]  

24. Y. Liang, H. Liu, F. Wang, H. Meng, J. Guo, J. Li, and Z. Wei, “High-efficiency, near-diffraction limited, dielectric metasurface lenses based on crystalline titanium dioxide at visible wavelengths,” Nanomaterials 8(5), 288 (2018). [CrossRef]  

25. M. Meem, S. Banerji, C. Pies, T. Oberbiermann, A. Majumder, B. Sensale-Rodriguez, and R. Menon, “Large-Area, high-numerical-Aperture multi-level diffractive lens via inverse design,” Optica 7(3), 252–253 (2020). [CrossRef]  

26. P. Lalanne and P. Chavel, “Large-area, high-numerical-aperture multi-level diffractive lens via inverse design: comment,” Optica 8(7), 1009–1010 (2021). [CrossRef]  

27. P. Wang, N. Mohammad, and R. Menon, “Chromatic-aberration-corrected diffractive lenses for ultra-broadband focusing,” Sci. Rep. 6(1), 21545–7 (2016). [CrossRef]  

28. S. Banerji, M. Meem, A. Majumder, B. S. Rodriguez, and R. Menon, “Imaging across an unlimited bandwidth: Is it possible?” in Optics InfoBase Conference Papers (2020).

29. N. Mohammad, M. Meem, B. Shen, P. Wang, and R. Menon, “Broadband imaging with one planar diffractive lens,” Sci. Rep. 8(1), 2799 (2018). [CrossRef]  

30. S. Zhang, A. Soibel, S. A. Keo, D. Wilson, S. B. Rafol, D. Z. Ting, A. She, S. D. Gunapala, and F. Capasso, “Solid-immersion metalenses for infrared focal plane arrays,” Appl. Phys. Lett. 113(11), 111104 (2018). [CrossRef]  

31. F. Presutti and F. Monticone, “Focusing on bandwidth: achromatic metalens limits,” Optica 7(6), 624–631 (2020). [CrossRef]  

32. J. Engelberg and U. Levy, “Achromatic flat lens performance limits,” Optica 8(6), 834–845 (2021). [CrossRef]  

33. S. Banerji and B. Sensale-Rodriguez, “A Computational Design Framework for Efficient, Fabrication Error-Tolerant, Planar THz Diffractive Optical Elements,” Sci. Rep. 9(1), 5801 (2019). [CrossRef]  

34. P. J. Van Laarhoven and E. H. Aarts, “Simulated annealing,” in Simulated annealing: Theory and applications (Springer, 1987), pp. 7-15.

35. J. Engelberg and U. Levy, “Standardizing flat lens characterization,” Nat. Photonics 16(3), 171–173 (2022). [CrossRef]  

36. E. Wolf and E. W. Marchand, “Comparison of the Kirchhoff and the Rayleigh–Sommerfeld Theories of Diffraction at an Aperture,” J. Opt. Soc. Am. 54(5), 587–594 (1964). [CrossRef]  

37. D. Faklis and G. M. Morris, “Spectral properties of multiorder diffractive lenses,” Appl. Opt. 34(14), 2462–2468 (1995). [CrossRef]  

38. X. Liu, X. Kong, and C. W. Qiu, “matlab_code_RPSA,” figshare (2022), https://doi.org/10.6084/m9.figshare.19597276