Impact of fabrication errors and refractive index on multilevel diffractive lens performance

Multilevel diffractive lenses (MDLs) have emerged as an alternative to both conventional diffractive optical elements (DOEs) and metalenses for applications ranging from imaging to holographic and immersive displays. Recent work has shown that by harnessing structural parametric optimization of DOEs, one can design MDLs to enable multiple functionalities like achromaticity, depth of focus, wide-angle imaging, etc. with great ease in fabrication. Therefore, it becomes critical to understand how fabrication errors still do affect the performance of MDLs and numerically evaluate the trade-off between efficiency and initial parameter selection, right at the onset of designing an MDL, i.e., even before putting it into fabrication. Here, we perform a statistical simulation-based study on MDLs (primarily operating in the THz regime) to analyse the impact of various fabrication imperfections (single and multiple) on the final structure as a function of the number of ring height levels. Furthermore, we also evaluate the performance of these same MDLs with the change in the refractive index of the constitutive material. We use focusing efficiency as the evaluation criterion in our numerical analysis; since it is the most fundamental property that can be used to compare and assess the performance of lenses (and MDLs) in general designed for any application with any specific functionality.

impact of fabrication errors and refractive index on multilevel diffractive lens performance Sourangsu Banerji, Jacqueline cooke & Berardi Sensale-Rodriguez * Multilevel diffractive lenses (MDLs) have emerged as an alternative to both conventional diffractive optical elements (DOEs) and metalenses for applications ranging from imaging to holographic and immersive displays. Recent work has shown that by harnessing structural parametric optimization of DOEs, one can design MDLs to enable multiple functionalities like achromaticity, depth of focus, wide-angle imaging, etc. with great ease in fabrication. Therefore, it becomes critical to understand how fabrication errors still do affect the performance of MDLs and numerically evaluate the trade-off between efficiency and initial parameter selection, right at the onset of designing an MDL, i.e., even before putting it into fabrication. Here, we perform a statistical simulation-based study on MDLs (primarily operating in the THz regime) to analyse the impact of various fabrication imperfections (single and multiple) on the final structure as a function of the number of ring height levels. Furthermore, we also evaluate the performance of these same MDLs with the change in the refractive index of the constitutive material. We use focusing efficiency as the evaluation criterion in our numerical analysis; since it is the most fundamental property that can be used to compare and assess the performance of lenses (and MDLs) in general designed for any application with any specific functionality.
There has been a significant interest in the scientific community in recent times to reduce the thickness and weight of lenses to enable miniature and compact optical systems 1 . Conventional refractive lens, which harness refraction to guide light, fails to satisfy both of this criterion, as it tends to be heavy and bulky owing to its increasing curvature with the increase in numerical aperture 2 . Unlike their refractive analogue, lenses, which exploit diffraction to guide incident light, have already been shown to be thin (thickness ~ 2λ) and lightweight. This ability to maintain a constant thickness is simply achieved by decreasing the local period of the diffractive optic 3 . To ensure constructive interference, each incident ray must now locally phase shift to compensate for the variation in its total optical path length to the focal plane. For conventional diffractive lenses, this is achieved by engineering the path traversed by the ray within the diffractive lens itself 3,4 .
In terms of performance, with respect to its refractive counterparts, traditional blazed, or diffractive lenses with almost optimal continuous phase distribution have already been shown to achieve 100% efficiency 5 . In principle, however, the analysis of such blazed gratings or diffractive lenses relies on the thin element approximation (TEA) which assumes no diffraction occur within the structure, i.e., it assumes that the grating will impose its structural shape on the incident beam wave-front. This is clearly not the case at high NA as diffraction occurs within the gratings with power being diverted to guided-mode resonances instead of the propagating modes 6 . In addition to this, conventional diffractive lenses also have poor broadband performance due to significant chromatic aberrations. The first problem is mitigated with parametric optimization of constituent elements of the diffractive lens, and the second problem is avoided to some extent through harmonic phase shifts 7 and by using higher orders of diffraction 8 . However, a harmonic diffractive lens (or multi-order diffractive lens) is, in principle, a hybrid refractive-diffractive lens. Such a refractive-diffractive based approach is limited only to a discrete number of operating wavelengths. Alternate methods to tackle the problems mentioned above involved harnessing appropriately designed pure phase masks to mimic the functionality of an achromatic lens or utilizing an axilens to achieve broadband achromatic focusing 9,10 . From such a perspective, multilevel diffractive lenses (MDLs), as shown in Fig. 1a 32,33 . We hypothesize that the achievable bandwidth with MDLs is "unlimited" for practical purposes and is only constrained by the quantum efficiency of the image sensor. Furthermore, we also showcased MDLs with a Field of View (FOV) up to 50˚ for wide-angle imaging 26 as well as MDLs with a Depth of Focus (DOF) imaging of up to 6 m in the NIR 34 . Apart from this, MDLs have also been utilized to create broadband holograms enabling multi-plane image projection 35 and in holographic displays for AR/VR applications 36 . Computational imaging with single and multi-aperture MDLs have also been showcased 37-39 along with its potential for applications in photovoltaics 40 .
Because of this stupendous progress made using MDLs, a common concern amongst most of the work in this field is the fact that there is still a discrepancy among simulated and experimental efficiency values. This ultimately affects the system performance. We acknowledge the fact that researchers indeed have made efforts to justify these discrepancies by resorting to fabrication errors, as is associated with a non-industrial grade fabrication facility. Yet still, we firmly believe that it would be of immense help to designers to understand how (and in which cases) fabrication errors do affect the performance of MDLs at the initial stage (i.e., even before putting them into fabrication) and numerically evaluate the trade-off between efficiency and initial parameter selection. This will not only help to make a judicious decision while choosing the initial parameters when designing a MDL, keeping in mind the fabrication capability at the corresponding location, but will enable one to predict (at least provide the ballpark estimates of) efficiency (with a certain degree of accuracy) even if a certain MDL is designed and tested out.
To facilitate with this process of providing a suitable metric, in this work, we seek to perform a statistical simulation-based study on MDLs (primarily operating in the THz regime) to analyse the impact of various imperfections (single and multiple) on the final structure as a function of the number of ring height levels. Figure 1b,c depict the two major cases, i.e., error in ring height and ring width, which is often commonly encountered in the fabrication of an MDL. Both these cases have been studied individually. Furthermore, an amalgamation of both errors presented simultaneously in the final MDL structure is analysed, to paint a more practical picture for narrowband and broadband operation. Furthermore, we also evaluate the performance of the same MDLs with variations in the refractive index. An aspect from which error (primarily at the modelling stage) can also occur (e.g. due to variations in density when 3D printing) and propagate during the optimization stage. It is essential to study the variation of performance with respect to material dispersion because most often, it is difficult to ascertain the exact optical constants (refractive index and absorption coefficient) of the material that is going to be used in the manufacturing process. In that case, MDL designers are left to rely upon values of materials already provided in the pre-existing literature while doing the design. In addition to this, many groups or fabrication facilities use modified hybrids of the same standard material, to which other groups do not have access too. A study pertaining to this specific problem is done in 41 pertaining to 3D printable materials used in manufacturing THz diffractive optical elements at a common fabrication facility versus the author's own laboratory. Finally, even if the design is done with a specific material; to have an idea of how well the same design may perform with a different material with similar properties will help one to gauge the robustness of the designed MDL structure.
The rationale behind the choice of focusing efficiency as the evaluation metric stems from the fact that focusing efficiency is the most basic criterion, which guides the field of lens design for any desired application with any tailored functionality. Other metrics like EDOF, FOV magnification, aberrations can also be included; but these are often not the principal criterion that one looks at when designing an optical element. The choice of performing this study on MDLs primarily operating the THz regime is due to two main reasons [42][43][44][45] . First, the ease in modelling larger unit cells ("rings") leads to the total number of elements in the entire structure to be www.nature.com/scientificreports/ small and tractable; hence, it is easy to capture the trade-off in efficiency with faster simulations and develop accurate prediction metrics. Therefore, the THz regime is an ideal frequency band to capture the real significance of this work 44,45 . Second, since the THz band lies in almost the centre of the electromagnetic spectrum, hence the metrics developed in this paper can be scaled up or down to both the microwave as well as optical ranges without any loss of generality. One can very well argue here that since the current study is specifically dealing with MDLs operating in the THz regime, which could be fabricated by 3D printing, that at visible and infrared wavelengths, where the MDLs are fabricated using other approaches, e.g., grayscale lithography or multi-step lithography, other error sources like misalignment and feature rounding 5 apart from the ones mentioned here will be significant. These are effects that indeed should also be analyzed and that are dependent on the particularities of the MDL fabrication processes. However, as shown, e.g., in [17][18][19]27,28,32,33 , the error in MDL's ring width and height, as analyzed in this study, still can contribute as a significant source of error to a great extent towards the reduction in MDL efficiency at these wavelengths. At optical wavelengths, the errors in height and width herein analyzed are referred to as "misetch" and "feature size error, " respectively 5 .
Finally, we would also like to highlight some limitations of the current study. First, the work, as reported here, only applies to MDLs operating under low NA (NA ≤ 0.3). At this point, we are not truly aware of how the error sensitivity varies for Mid-(NA > 0.3 and NA ≤ 0.6) to High-NA (NA > 0.6) MDLs. This is because for Mid-to High-NA MDLs, due to "shadowing effect, " one would be required to perform full-wave simulations instead of simply relying on scalar diffraction theory for accurate estimates. Second, the analysis (to be specific, the scalar model) in this study, neglects reflection losses at the interface between PLA and air. This, however, is not expected to significantly alter the overall outcome of the analysis, as discussed in previous works, e.g. 27,28,32,33 . Third, the study is limited to the sources of error already mentioned, i.e., error in ring width and ring height. In order to provide even more accurate predictions, one would be typically required to incorporate additional error and modeling constraints on top of this existing framework. To re-iterate the same statement made earlier, in this work, our focus is to develop a simple metric which will enable one to predict (at least provide the ballpark estimates of) efficiency (with a certain degree of accuracy).

Results and discussion
The rotationally symmetric MDLs were designed via a non-linear search method, namely Gradient Descent Assisted Binary Search (GDABS) algorithm coupled with scalar diffraction theory under Fresnel approximation. Full description of the design process is already explained in 29,30 . Therefore, we choose to omit an in-depth discussion of the same here. However, to surmise, each of the designed MDL consists of concentric rings of width equal to a pre-defined value with the ultimate objective of maximizing the focusing efficiency across its bandwidth of operation. Radial symmetry of these structures is exploited to speed up the computation and reduce the optimization time. For our current study, we designed five different MDLs for both narrowband (0.2 THz) as well as broadband (0.1 THz to 0.3 THz) operation. Therefore, in total, ten different MDLs were designed. The ring height distribution and the relevant geometric design parameters are provided in the Supplementary Information.
Speaking of geometric design parameters, the design values were chosen, keeping in mind the fabrication constraints associated with any off-the-shelf hobbyist 3D printer available in the market. To be specific, the parameters were chosen, keeping in mind the lowest vertical and lateral resolution possible with the 3D printers available from Cura Ultimaker, Dremel DigiLab, MakerBot Replicator, and Formlabs range of 3D printers. Knowledge of vertical and lateral resolution possible with this range of 3D printers is of prime importance because vertical resolution dictates the amount of error in pixel height, whereas lateral resolution dictates the error in pixel width of the fabricated MDL. Each MDL was designed to be 24 mm in diameter. The MDLs designed to operate at a single frequency of 0.2 THz comprised of multilevel concentric rings having a maximum thickness (h max ) = 1.6 mm whereas the broadband MDLs have a maximum thickness (h max ) = 3.2 mm. Each ring of the designed MDLs has a width (w) = 0.4 mm. Therefore, the total number of such rings within each MDL = 30. For the purpose of this study, the number of distinct ring height levels were varied from P = [8,16,32,64,128] which dictated the minimum thickness (h min ) = [0.2 mm, 0.1 mm, 0.05 mm, 0.025 mm, 0.0125 mm] for narrowband MDLs and (h min ) = [0.4 mm, 0.2 mm, 0.1 mm, 0.05 mm, 0.025 mm] for broadband MDLs. The focal length was fixed at f = 50 mm, which translates to a numerical aperture (NA) of 0.2334. The material chosen for the design of the MDLs was PLA, and its dispersion values (which was inputted into the optimization algorithm) are provided in the Supplementary Information. The reason behind choosing PLA is due to its widespread availability as a common 3D printable material and its negligible loss within the frequency of operation of the designed MDLs. However, other materials could also have been used.
A statistical standard deviation-based simulation error model like in 17,32,33 was undertaken to study the impact of fabrication error (ring height and ring width) due to the following reasons: One, the standard deviation is always considered in relation to the mean (or average) since the mean by itself is usually not very useful. Two, the standard deviation is the best measure of variation since it is based on every item in the distribution 46 . Third, the standard deviation is less affected by fluctuations of sampling than most other measures of dispersion. Fourth, the standard deviation is most prominently used in carrying out further statistical studies like computing skewness, correlation, etc. 47 . Therefore, the use of a standard deviation-based error model seems justified. In this study, we varied the ring height (Δherr) as well as the ring width (Δwerr) of the optimized MDLs first individually and then as an amalgamation of both and re-calculated the focusing efficiency. The error terms for ring heights and ring widths were embedded in the pupil function T x, y, of the MDL which can be written as: Similarly, one can derive the same expression for the Δwerr(m), whose expression can be written as, Once the error terms were incorporated into the pupil function, the beam propagation from the lens plane to the focal plane was carried out using the Fresnel-Kirchoff diffraction integral shown below: where f is the focal length, k equals 2π , (x, y) are coordinates in the lens plane, ( x ′ , y ′ ) are coordinates in the focal plane and T x, y, is the pupil function of the MDL. We assume that only the light that falls inside the MDL's active area is diffracted, and whatever falls outside this area propagates through unaltered. The intensity in the focal plane is given by: Finally, the focusing efficiency can be expressed as: where m equals the number of wavelength samples, variables ( w xj , w yj ) are the full width half maximum (FWHM) of the theoretical diffraction-limited PSF in both the x-and y-direction for the respective wavelength sample. In the numerator, I f x ′ , y ′ , j is the intensity at the focal plane; whereas, in the denominator I i x, y, j represents the intensity at the lens plane (rather total power impinging on the MDL). The fixed quantity 3 2 in the term 3w xj 2 denotes that the optimization routine will try to maximize the intensity within a spot of diameter equal to 3 times the theoretical diffraction-limited FWHM of the PSF at the focal plane. For each number of distinct ring height value P, we recorded ten sets of observation data. From this observation data, we then proceeded to calculate the mean as well as the 90% confidence interval around this mean. We then plotted the same in Figs. 2 and 3. Figure 2a,b depicts the impact on average focusing efficiency for the first fabrication error scenario, i.e., standard deviation-based error in ring height for the designed MDLs under narrowband (0.2 THz) and broadband operation (0.1 THz-0.3 THz) respectively. We computed the focusing efficiency of the MDLs as the power within a spot of diameter equal to 3 times the FWHM of the spot divided by the total power incident on the lens 3 . Later we averaged the individual focusing efficiencies over the total number of design frequencies. This is exactly what Eq. (8) is in mathematical notation. On similar lines, Fig. 2c,d portrays the impact on average focusing efficiency second fabrication error scenario, i.e., standard deviation-based error in ring width for the designed MDLs under both narrowband (0.2 THz) and broadband operation (0.1 THz to 0.3 THz) correspondingly. We would like to remind the readers that the terms "average focusing efficiency" and "average efficiency" are the same and will be used interchangeably in this discussion. Moreover, the term "average efficiency" for the narrowband MDL is just "efficiency" since the number of frequency samples in this case = 1. Moving ahead, we considered the standard deviation in the error of both ring height and width up to 500 µm because the fall in efficiency for all the plots is ~ 50%, barring only the first case for narrowband MDLs operating at a single frequency.
Some general observations from the plots in Fig. 2 are as follows: The MDLs designed to operate at only a single frequency are very resilient to errors in both ring height and width as compared to its broadband counterparts. This is easy to understand given the fact that the constraint of focusing a larger number of frequencies at a single point in the observation plane requires a higher degree of control by the diffractive structures, and hence a slight mismatch would significantly degrade its performance. An important corollary to the former (2) T x, y, = 1 + m circ(r − m(� w + � werr )) e ia( )(� h +� herr )pm − 1 .  [8,16,32,64,128]; whereas the width of each ring is kept constant at 0.4 mm during the entire optimization. Even if one assumes that the optimized MDL might converge to a solution where the height of two adjacent rings have the same value such that both the rings can now be considered as a single ring of 0.8 mm thickness, the differential width (Δw) still has a smaller degree of freedom in contrast to its height. A smaller degree of freedom here refers to the fact that the number of discrete width values which the ring width can assume during the optimization is smaller than the number of discrete height values which the ring height can assume. Therefore, as evidenced from Fig. 2, this ultimately leads to an optimized MDL structure, which is much more resilient to errors in height than in width. Finally, to quantify the results, for a 20% reduction in focusing efficiency, the narrowband MDLs require a ~ 350 µm error in ring height as compared to a ring width error of ~ 250 µm. For the broadband MDLs, the values are ~ 250 µm and ~ 200 µm in ring height and width, respectively. The "20% decrease in focusing efficiency" metric has been adopted from 3 , where it is seen that the change in the PSF plots is prominent for this value of decease in focusing efficiency.
A second set of ten MDL designs were also undertaken to understand the influences of fabrication error on other geometries generated by the algorithm (with the same design parameters as those discussed herein). The MDL pixel profiles and the corresponding error analysis plots, following the analysis in Figs  www.nature.com/scientificreports/ gauge how the point spread functions (PSFs) deteriorate as errors are augmented. This will help one to quantify the basic imaging quality of an optical system. We also analysed the relative error arising from the inner (central) and outer (periphery) rings of the MDL to investigate how sensitive is the performance of the structure to localized variations taking place only in certain regions of the MDL. For this, we choose to perform the same standard deviation based error in pixel heights analysis and its impact on the focussing efficiency, but now we divided the total number of rings (i.e., 30) within the MDL into two equal sections, i.e., inner and outer regions consisting of 15 rings each. The plots for the same are shown in Fig. S8 of the Supplementary Information. We performed this analysis on MDLs designed for both the narrowband ( Fig. S8 (a,b)) as well as the broadband (Fig. S8 (c,d)) operation. There are two main key takeaways. First, the impact of error on efficiency is higher for the outer rings as compared to the inner rings. A qualitative 5%-10% relative error is observed across the plots for both narrowband as well as the broadband operation. This is in line as to what one can expect, for instance, for any DOE, where it becomes difficult to focus (diffract) light from the outer zones (due to larger spreading) in contrast to the central zones 48,49 . A more rigorous treatment of the same is provided in 50 . Second, the impact on focussing efficiency was higher for the broadband MDLs (~ 7-10%) in comparison to the narrowband MDLs (3-6%). This could be understood from the perspective that the challenge to achieve achromatic broadband focusing is much higher than focusing single frequencies.
Finally, we study the impact of change in the refractive index of the constitutive material on the average focusing efficiency of the MDL structures, again, at both a single frequency of 0.2 THz as well as a broadband regime from 0.1 THz to 0.3 THz. The respective plots for both these cases are depicted in Fig. 4a,b, respectively. A point of note here is that the absorption coefficient is very negligible in this frequency range, i.e., 0.1 THz to 0.3 THz, as is evidenced from the plot in Fig. S2 of the Supplementary Information. When the refractive index values were changed during this study, the term a(λ) in the pupil function T x, y, of the MDL in Eq. (2) was updated accordingly. The term a(λ) can be written as k[{(n + �n) + iK } − 1] where n ∈ [−0.2, 0.2].
From Fig. 4b it is observed that the fall in the average focusing efficiency is steeper for the broadband MDLs in contrast to the narrowband MDLs in Fig. 4b. The narrowband MDLs are also more resilient to the change in refractive index too. This, again, can be attributed to the "potential challenge" in handling more design constraints with accuracy. Another important observation that is critical here is the fact that for the change in refractive index ( i.e. n ∈ [−0.2, 0.2] ), the MDL solutions for both narrowband and broadband become equally unstable. This tells us that irrespective of the bandwidth or the number of ring height levels, the optimized structures will just not work if the dispersion values are way off than what was used during the design stage. In conclusion, we have performed a statistical simulation-based study on MDLs designed to operate in the THz regime to analyse the impact of various fabrication imperfections in the final structure as a function of the number of ring height levels. In addition to this, we also analysed the performance of these same MDLs with the change in the refractive index. We firmly believe that this work provides important information to MDL designers and researchers alike to compare and assess the performance of MDLs as well as offer fundamental insight into designing highly efficient and robust MDLs with given fabrication constraint.