Hexagonal diffractive optical elements

: Diffractive optical elements (DOEs) have widespread applications in optics, ranging from point spread function engineering to holographic display. Conventionally, DOE design relies on Cartesian simulation grids, resulting in square features in the final design. Unfortunately, Cartesian grids provide an anisotropic sampling of the plane, and the resulting square features can be challenging to fabricate with high fidelity using methods such as photolithography. To address these limitations, we explore the use of hexagonal grids as a new grid structure for DOE design and fabrication. In this study, we demonstrate wave propagation simulation using an efficient hexagonal coordinate system and compare simulation accuracy with the standard Cartesian sampling scheme. Additionally, we have implemented algorithms for the inverse DOE design. The resulting hexagonal DOEs, encoded with wavefront information for holograms, are fabricated and experimentally compared to their Cartesian counterparts. Our findings indicate that employing hexagonal grids enhances holographic imaging quality. The exploration of new grid structures holds significant potential for advancing optical technology across various domains, including imaging, microscopy, photography, lighting, and virtual reality.

Unfortunately, Cartesian grids have several practical disadvantages.First, they exhibit a high degree of anisotropy in the plane, with diagonal frequencies represented more poorly than axially aligned frequencies [12,13].For this reason, deformable mirrors and other phase modulators intended for adaptive optics often opt for a hexagonal packing of the actuators.Due to its six-fold symmetry, this hexagonal arrangement exhibits better isotropy [14].
In recent years, hexagonal pixel arrangements have shown potential for achieving better resolution and efficiency in fields such as image processing and computer graphics [14][15][16].In the field of sensing techniques, hexagonal sampling gives wider spectra than rectangular with the same pixel density [12].Thanks to its six-fold symmetry, the hexagonal lattice offers more isotropy in representing complex elements than Cartesian grids.
Another disadvantage of using Cartesian grids for DOE design is that they yield square-shaped features, placing a significant burden on the accuracy of the fabrication process.For instance, when using photolithography, it is common to experience some corner rounding on the features [17], leading to deviations between the fabricated DOE and its intended design.Moreover, in Cartesian grids, features may touch on a corner, which results in topological changes during fabrication, particularly when using equipment with limited resolution capabilities.As an example, Fig. 1 shows the pattern produced by a mask writer (Heidelberg DWL2000), exhibiting corner rounding and topology artifacts.In contrast, hexagonal features closely approximate the diffraction-limited spot of common lithography equipment such as mask writers, and adjacent features in a hexagonal grid always share an edge, which preserves the topology of the design.These advantages potentially allow us to design for smaller feature sizes, and corresponding stronger diffraction or higher diffraction efficiency using the same lithography equipment.In optics, the use of hexagon arrays has been studied in the field of integral imaging, lens arrays [17][18][19], image display [20], sensors [21], and phase modulators [22,23].In this work, we explore the utilization of hexagonal grids for DOE design, which combines the rapid Fourier-based simulation of light propagation from Cartesian grids with the sampling and fabrication advantages of hexagonal grids.This combination is made possible through the use of a hexagonal Fast Fourier Transform (HFFT), which can be considered a drop-in replacement for the Fast Fourier Transform (FFT) in any existing DOE optimization framework, be it GS or a more recent gradient back-propagation method.While our approach is applicable to various DOE design tasks, our primary focus lies in holography applications.Holography offers a straightforward means of characterizing the quality disparities, both in terms of simulation accuracy and direct measurements of fabricated DOEs.

μm 5 μm
Our contributions encompass several key aspects.Firstly, we adopt a novel 2D hexagonal arrangement for DOE geometries, opting to utilize the Array Set Addressing (ASA) coordinate system [15] for data storage and processing.This innovative approach forms the basis for our subsequent work.Secondly, we implement hexagonal scalar diffraction with differentiable functions, facilitating the back-propagation of gradients.This capability enables us to perform hexagonal Fresnel diffraction simulations and conduct inverse design.To validate our methodology, we compare the simulation accuracy between hexagonal and rectangular sampling schemes.Furthermore, we apply our techniques practically by converting images, serving as holograms, into hexagonal-sampled images.We utilize back-propagation and GS to inversely derive their wavefront, demonstrating the real-world applicability of our approach.The DOEs are fabricated using a combination of photolithography and reactive ion-etching (RIE) techniques, and subjected to comprehensive testing, with comparisons made against DOEs featuring Cartesian grids.As a result, practical holographic examples empirically validate the correctness and effectiveness of our implementation.To facilitate further research in the community, we will release the associated code online at https://github.com/vccimaging/HexDOE.

Methods
To demonstrate the benefits of the hexagonal grid in DOE design, we first introduce general data processing operations and fast Fourier transforms for the hexagonal grids in the ASA format following [15].Next, we illustrate how to implement the Fresnel diffraction method in the hexagonal grids for general wave propagation.Last, for the purpose of demonstration, we implement the classic GS algorithm [9] as well as a differentiable back-propagation method for the inverse design of a holographic hexagonal DOE.

Hexagonal coordinate system and HFFTs
Among multiple choices for hexagonal grid representation [24][25][26][27][28], the ASA coordinate system [15] is best suited for data storage and allows easy implementation of algorithms.As shown in Fig. 2, two rectangular grids are employed to represent a hexagonal grid.The hexagons are identified using three coordinates, where the coordinates represent the array (a), row (r), and column (c), respectively.
Fig. 2. Hexagonal coordinate system: Array Set Addressing.The data on hexagonal grids are stored as two interleaved rectangular arrays.Each hexagon pixel in the 2D space can be uniquely represented by (a, r, c) coordinates, where a is either 0 or 1 to address which rectangular array the pixel is stored in.r and c address the row and column respectively.s h is the horizontal spacing of the hexagonal grids, and e is the edge length of a single hexagon.
Note that we adopt the pointy-top-orientated [29] regular hexagon layout in the ASA coordinate.The horizontal spacing between neighboring hexagons s h measures the sampling density of the hexagonal grids, and it follows that s h = √ 3e, where e is the edge of the hexagon.The forward and inverse hexagonal discrete Fourier transforms in ASA coordinates [30] can be expressed as where f (a, r, c) is the spatial signal, F(b, s, d) is its Fourier transform, and Throughout the equations, both spatial and frequency dimensions are 2 × n × m (n is vertical and m is horizontal), and a, b Notice that Eq. ( 2) is separable and can be further expressed in the form of 1D fast Fourier transforms, where ]︃ . ( These expressions can be realized using standard 1D fast Fourier transforms.We implement the HFFT using PyTorch [31] to support hardware acceleration and automatic differentiation. Similar to the 2D Cartesian data, the HFFT assumes that the signal is periodic.To perform the Fourier transform properly, it is necessary to apply a Fourier shift before FFT.We exploit the periodic property of the hexagonal data as shown in Fig. 3(a) to derive proper shift operation.The illustrated hexagonal Fourier shift is used to align the 2D signal correctly.The periodic hexagon sampling domain can be shifted to an interlocking rectangular periodic domain.
To ensure isotropy in the frequency domain, we let 3n = 2m.Otherwise, the hexagonal grid in the frequency domain will be scaled unequally along vertical and horizontal directions.This also gives the spatial sampling region to be regular hexagon.
In Fig. 3(b), we show the whole hexagonal Fourier transform process of a square signal that uses the HFFT.In the spatial domain, we first sample and shift the data according to the (a, r, c) coordinates.The HFFT is then applied to compute the discrete Fourier transform, followed by an inverse hexagonal Fourier shift to restore to the original order.

Hexagonal scalar diffraction
Wave propagation is commonly expressed by the Fresnel diffraction [11] in the Cartesian coordinates.Similarly, we can employ the forward and inverse HFFTs to rewrite the Fresnel diffraction in the hexagonal coordinates, where F z {•} is the Fresnel propagation, z is the propagation distance, r is a 2D vector that represents the position at the source and observation plane, u 0 (r) and u z (r) are complex fields at the source plane and image plane respectively, and h z (r) is the impulse response given by h z (r) = e jkz jkz exp where k = 2π/λ is the wave number, and |r| is the distance to the z-axis at the observation plane.
For the hexagonally-sampled field with the ASA representation and a spacing of s h , we have The scalar diffraction for a DOE can be mathematically modeled as a complex linear system.This system takes a plane wave as input and generates the diffraction field at a given propagation distance z as output.The DOE phase diagram disrupts the input wave by adding a phase shift to it.Our objective is to extract the phase information from the output intensity of the system.This can be accomplished by minimizing the squared error between the reconstructed and target intensities.

Inverse hexagonal DOE design
The process of inverse design of the DOE is illustrated in Fig. 4. The first step in this process is to retrieve the phase distribution from the desired image by reversing the forward simulation of Fresnel diffraction.This phase distribution is then used to create the wavefront of the DOE with some fabrication constraints.With this, we can design the DOE to produce the desired hologram.
The problem we are trying to solve can be expressed as follows: where ϕ is the phase of the DOE to be optimized, u 0 (r) is the incident plane wave, and I target is the intensity of the hologram we wish to reconstruct.
With the basic arithmetic operations and Fourier transforms for the ASA system, any algorithms developed for Cartesian coordinates can readily be modified as their hexagonal counterparts to solve the DOE optimization problem.In this work, we implement two algorithms, GS and back-propagation, for the purpose of demonstration.Other methods can be implemented similarly.
The GS algorithm is an iterative phase retrieval algorithm composed of forward and backward calculations.With the hexagonal Fourier transform and its inverse, we are able to realize the update of phase by replacing the amplitude of backward calculation with the target amplitude.
The other algorithm is back-propagation.The inverse problem in Eq. ( 9) can be efficiently solved using advanced iterative algorithms, such as stochastic gradient descent [32] or Adam [33].The optimization algorithm is straightforward and can be implemented using state-of-the-art deep learning modules.
Compared to traditional iterative methods such as GS, using back-propagation is more flexible.The advanced optimization takes advantage of avoiding the local minimum [32] to increase the accuracy and provides the use of regularization terms [34] to add constraints to the solution.Regarding the phase retrieval problem, various results with different requests could be obtained such as smoothed intensities, concentrated energy distribution, adjustable aperture size, etc.

Simulation
We consider a holography system with a hexagonal DOE.The target image size is 4 mm × 4 mm and the wavelength of the source plane wave is 520 nm.Zero padding is added around the DOE to ensure the simulation accuracy [35], which gives the whole simulation area to be 8 mm × 8 mm.We use the same configuration for both Cartesian and hexagonal grids for a head-to-head comparison.The horizontal and vertical Cartesian sampling spacings are both s r .To ensure identical sampling density in the hexagonal coordinates, we must let A series of rectangular sampling spaces ranging from 1 µm to 10 µm have been used to test the fabrication limitation.For the sake of convenience, we use the term "feature size" ∆r, which is equal to the rectangular sample spacing s r , to represent the hexagonal and rectangular resolution with the same sampling density.For example, for the feature size of 1 µm, the rectangular sampling spacing s r = 1 µm, and the hexagonal sampling spacing s h ≈ 1.075 µm.

Forward simulation accuracy
To measure the simulation accuracy, we first compare the diffraction results of a square aperture to the analytical solution of Fresnel diffraction [11].Various widths for the square aperture are used and the propagation distance is 70 mm.To test the sampling accuracy for tilted features, a slightly rotated square beam (7°) is also used for comparison with the one without rotation.The feature size is 5 µm and the width of the aperture ranges from 0.1 mm to 2.0 mm.
The results are illustrated in Fig. 5.The peak signal-to-noise ratios (PNSRs) of the simulations are compared in Figs.5(a) and 5(d).For an aperture of 1.0 mm width, the corresponding simulation errors are shown in Figs.5(b) and 5(c) for the hexagonal and Cartesian grids, respectively.It is clear that hexagonal sampling can reduce the error along the x-axis due to its interlaced geometry.However, even though hexagonal coordinates exhibit superior accuracy in these two scenarios, the results diverge when the rotation angle surpasses approximately 30°as in Fig. 5(d).The accuracy of Cartesian coordinates fluctuates significantly with varying rotation angles.This observation underscores that, for this square aperture, Cartesian coordinates are highly sensitive to their direction, whereas hexagonal coordinates demonstrate much greater robustness.A more equitable comparison with the circular beam is presented in Fig. 5(e), where accuracy remains uninfluenced by the sampling direction.In this context, hexagonal coordinates outperform Cartesian coordinates across most radius values.It is also worth noticing that oscillation occurs when the aperture size varies due to the misalignment between the aperture and the sample grids.

Inverse design
To evaluate the performance, we choose various target images and inverse design the DOEs to generate corresponding holograms.For back-propagation, we choose the Adam optimizer with a starting learning rate of 0.01 to update the phase for each iteration.We start with a uniform distribution to achieve a smoother phase that avoids high-resolution features that hurt the fabrication result.The optimization ends after 100 iterations.
In Fig. 6(a), we show three examples, the KAUST logo, a photograph of Vermeer's "Girl with a Pearl Earring", and a spiral pattern (from top to bottom).The corresponding simulated intensity images, with ∆r = 10 µm, are shown in Fig. 6(b) for visual comparison.
In Fig. 7, we compare the imaging resolution of hexagonal and Cartesian designs using the USAF 1951 resolution chart.It contains features of different sizes, which can challenge the sampling.We first use linear interpolation to convert the original images to hexagonal-sampled data.The source is set to be a plane wave with a square aperture of 4 mm × 4 mm, which is the same as the image size.The propagation distance is 200 mm.
The target resolution chart and a zoom-in region are shown in Fig. 7(a).The simulated intensity images for hexagonal designs are shown in Fig. 7(b) and Cartesian designs in Fig. 7(c outcomes for hexagonal and Cartesian designs exhibit a high degree of similarity, attributable to their identical sampling density.The specific nuances diverge due to the stochastic nature of the optimization processes and variations in their sampling approaches.For example, when considering ∆r = 10 µm, the hexagonal design demonstrates superior performance in Group 2, whereas with ∆r = 5 µm in Group 1, the results are reversed.As for ∆r = 1 µm, the sampled targets closely align since the features are significantly larger than the spacing, resulting in phase retrieval outcomes that are exceptionally close.

Experimental results
The DOEs are fabricated by a 4-step workflow of consecutive photolithography and RIE steps as employed in [1,2,8].This yields a 16-level quantization of the DOE profile to approximate the continuous design.A chromium aperture is added to restrict stray light from entering the sensor.As illustrated in Fig. 8, to measure the hologram created by the DOE, a laser (Thorlabs PL201) with a wavelength of 520 nm is used as the source.A collimator (HONC HC-550) is used to create the plane wave.An image sensor (Lucid TRIO54S-CC) with a resolution of 2880 × 1860 and pixel size of 3.0 µm is used to collect the diffracted pattern of the DOE.implementation of phase retrieval algorithms is correct.The phase retrieval result of the GS algorithm is also fabricated and measured in Fig. 9(c) as a baseline reference.
To further challenge the fabrication resolution and explore the advancement of hexagonal grids, we reduce the aperture as well as the propagation distance to increase the diffraction angle.The aperture is reduced to a circular one with 2 mm diameter.The propagation distance is reduced to 70 mm.The results are illustrated in Fig. 9(d).Different from that in Fig. 9(b), the error between simulation and measurement is significantly increased due to the increased design difficulty.
To experimentally compare the performance of hexagonal and Cartesian designs, we show the hexagonal design in Fig. 10  Note that in the middle of the spiral pattern, the hexagonal design shows more details than the rectangular design.To quantitatively compare the performance, we employ the PSNR as a metric.For hexagonal grids, PSNR = 6.155 and SSIM = 0.0728, and for rectangular grids, PSNR = 5.887 and SSIM = 0.0591.Although the distinction between hexagonal and Cartesian becomes subtle with a small diffraction angle, our findings indicate incremental improvements both in visualization and statistical metrics for this challenging scenario.

Discussion
We have successfully implemented the wave propagation simulation and inverse design algorithms on the hexagonal grids.The hexagonal sampling gives a better simulation accuracy due to its higher degree of isotropy.The measurement of the fabricated DOEs shows that the reconstruction quality using a hexagonal grid is superior to that of Cartesian grids.
Although we have demonstrated that the fundamental Cartesian calculations and optimization algorithms for diffraction simulation can be transferred to a hexagonal grid, there are still complex design strategies that remain explored.Since traditional 2D data is mostly processed in Cartesian coordinates, when converted to hexagonal coordinates, they often become more complicated.
In terms of computational efficiency, the 2D HFFT consists of two 1D FFT processes which takes more time than a direct 2D FFT for Cartesian grids.Moreover, HFFT requires additional operations that certainly consume more computational resources.This is shown in execution time.For example, for back-propagation of a hexagonal grid with a data size of 2×462×693 = 640, 332, it takes 31.325seconds to finish the optimization with 100 iterations on CPU whereas for a Cartesian grid with a size of 800 × 800 = 640, 000 it takes 7.825 seconds.To improve the calculation efficiency, a GPU implementation of HFFT could significantly enhance efficiency in the future.
Our experiments show the advantages of the hexagonal grid representation both in simulation and in the final fabricated results.The improvements in simulation can be attributed to the more isotropic sampling properties of the hexagonal grid.The additional benefits, while present, are subtle in the fabricated hologram results, demonstrating a nuanced improvement in matching the designed DOE and its fabricated counterpart with the hexagonal coordinate.It should be noted, however, that these additional benefits are modest, pointing to other sources of artifacts that affect the fabrication process, such as uneven etching rates in the RIE process, which equally affect both methods [36].Moreover, it's essential to highlight that the applications of DOEs go beyond holography, potentially capitalizing more on the advantages of hexagonal grids.For the existing optical elements such as hexagonal-packing SLMs and deformable mirrors that employ the hexagonal grids, the hexagonal wave propagation approach provides a more direct and efficient way to simulate and design, compared to resampling to Cartesian grids that introduces error.The hexagonal representation can be extended to advanced technologies, providing distinctive benefits in applications such as microscopic imaging and laser beam shaping, where sensitivity to fabrication intricacies is particularly pronounced.

Conclusion
In this work, we have introduced the design and optimization of DOEs in the hexagonal coordinates.The scalar diffraction formula in the form of hexagonal coordinates is derived and tested with high simulation accuracy.We show that the general algorithms used in conventional Cartesian designs can easily be transferred to hexagonal designs with the proposed methods.As a demonstration, we have implemented two phase retrieval algorithms, including GS and back-propagation, to recover the target holographic images with hexagonal sampling.DOEs with hexagonal grids are fabricated and measured to show the benefits of using hexagonal representation.Compared with traditional Cartesian grids, hexagonal grids reduce the edge effects thanks to the 6-fold symmetry.Due to higher uniformity along different directions, the hexagonal design has the potential to create more detailed diffracted patterns.We envision that hexagonal DOE offers an alternative design strategy for diffractive optics with applications that are not limited to holography, but also beam shaping, PSF engineering, and other advanced imaging and photography techniques.
Funding.King Abdullah University of Science and Technology (Individual Baseline Funding).

Fig. 1 .
Fig. 1.Fabricated masks with hexagonal and square patterns by laser direct writing (Heidelberg DW2000).Square features are more affected by corner rounding artifacts (areas shaded in green) and touching corners (areas shaded in red), while hexagonal features better preserve the geometry.Microscopic images are taken by a Nikon Eclipse L200N microscope.

Fig. 3 .
Fig. 3. Illustration of Hexagonal Fourier shift and Fourier transform.(a) Fourier shifts in the spatial domain and frequency domain.(b) Workflow of hexagonal Fourier transform by HFFT.

Fig. 4 .
Fig. 4. Inverse design of a hexagonal DOE.(a) The phase distribution of a target image is retrieved by iterative forward and backward propagation.(b) The optimized phase profile is quantized into 4 bits (16 levels) for fabrication.Finally, the fabricated DOE is used to produce the desired hologram.

Fig. 5 .
Fig. 5. Simulation accuracy comparison of Cartesian and Hexagonal grids with Fresnel diffraction.(a) PSNR of the simulated complex field compared to the analytical results using a square aperture.(b) Error map of the simulation result using a square beam with a width of 1.0 mm.(c) Error map of the simulation result using a rotated (7°counter-clockwise) square beam with a width of 1.0 mm.(d) PSNR of the rotated square beam with different rotation angles and a width of 1.0 mm.(e) PSNR of a circular beam with different aperture size
(a) and rectangular design in Fig. 10(b) for the spiral pattern with ∆r = 5 µm.

Fig. 10 .
Fig. 10.Comparison of the hexagonal and Cartesian measurements with ∆r = 5 µm.The same target image is used for hexagonal (a) and Cartesian (b) grids, respectively.Microscopic images are taken by a Nikon Eclipse L200N microscope.