## Abstract

The redistribution of an incoming radiation into several beams is necessary in telecommunication to demultiplex data signals. In the terahertz spectral range, it can be realized by easy-to-manufacture diffractive optical elements (DOEs) allowing to focus the radiation into multiple focal spots in a single plane. In this article, we present diffractive optical elements focusing THz radiation into three focal spots. Different focal spot distributions (symmetric and asymmetric) are designed using an iterative algorithm. The phase distribution forming asymmetric focal spots can be realized by iterative design, which is a novel approach, to our knowledge. Then, the structures are manufactured using a sintering-based 3D-printing method from polyamide 12 (PA 12) and measured in an experimental setup for 150 GHz frequency. A novel approach based on neural networks (NNs) is proposed to optimize the phase delay maps of the structures to further improve their performance – the higher efficiency and the lower unwanted background noise.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Application of diffractive optical elements (DOEs) for the terahertz spectral range benefits from the combination of the vast flexibility of beam shaping with the simplicity of manufacturing [1]. It is commonly known that the smallest details of the diffractive structures should be in the subwavelength range to assure their proper functioning. In the near-visible spectral range, it implies the usage of sophisticated fabrication techniques like laser or electron beam writing [2]. However, in the sub-THz region the wavelengths are in the order of at least tenths of a millimeter and much simpler methods, like additive manufacturing, can be used [3–5]. In this work, we focus on the methods of redirecting and focusing of the THz radiation into three spots in a single plane. The preliminary work has already been performed and discussed in [6]. It describes the optimization of the height and the shift of the diffraction grating to obtain only three diffraction lines ($-1^{st}$, $0^{th}$, $1^{st}$). When combined with the converging lens, such a composed structure allowed to obtain three focal spots at a chosen distance. Nevertheless, the higher orders of diffraction were still visible as the design was based on the binary grating. Another approach is the use of properly designed Dammann gratings [7,8] also combined with a converging lens. However, here the undesired focal spots in the designed matrix are usually observed [6]. In this work, we propose the application of iterative algorithms to mitigate higher orders of diffraction. Such algorithms offer one additional and non-trivial advantage which is the possibility of designing asymmetrical distributions of points. It is especially appealing when the increase of the bandwidth is of concern, like e.g in tomography systems [9,10]. In case of the previously described symmetrical 3-focal-point distributions (central $0^{th}$ order spot and two $\pm 1^{st}$ order spots, one shifted to the left and one to the right) the distances from both side focal spots to the center of the lens or the detector are equal. It means that signals from $\pm 1^{st}$ order focal spots reach the detector at the same time and are therefore indistinguishable. On the other hand, one can imagine an asymmetrical 3-point design, based on $0^{th}$, $1^{st}$ and $2^{nd}$ diffraction orders, all shifted to the same side in relation to the optical axis. The times of flight from all points to the detector will be then different and the proper demultiplexing of the signals will become possible. The application of DOEs realizing asymmetrical matrices of points in the THz spectral range is a novel idea discussed in this article. Additionally, optimization algorithms based on neural networks (NNs) are proposed for further improvement of the structures.

## 2. Design of the structures

In this article, we present a diffractive terahertz lens focusing the beam into three spots in a single focal plane. Two approaches have been investigated: symmetrical (two $\pm 1^{st}$ order focal points on both sides of the central one) and asymmetrical ($1^{st}$ and $2^{nd}$ order focal spots placed on one side of the central one). In the symmetrical approach, the designed structure will also be symmetrical, which facilitates future usage of the lens and ensures more even power distribution between particular focal spots as well. On the other hand, the asymmetrical approach allows to distinguish the signals reflected from or diffused by the sample by the time of arrival (different actual distances from particular focal points to the center of the lens or a detector). It opens up to a whole variety of possibilities for the multipoint scanning of the samples.

Diffractive structures and theoretical light distributions have been simulated with the modified convolution method of light propagation [11]. Additionally, an iterative algorithm, similar to Gerchberg-Saxton (G-S) [12] or ping-pong [13], has been used for the optimization of the structures. In-depth description of the optimization algorithm can be found in previous reports [14]. All simulations have been performed on 4096x4096 matrices of 117 µm pixels. This way, more than sufficient buffer region is provided around the structures and focal spots to avoid any side effects from the square-shaped calculation region. The pixel size is connected with the resolution of typically used 3D-printers [3,15]. The phase distribution of the designed structures and the calculated focal spot intensity distributions are illustrated in Fig. 1.

The presented DOEs have 200 mm of diameter, the focal length and the propagation distance have been set to 100 mm, while the frequency of the illuminating radiation is 150 GHz, which corresponds to a wavelength of 2 mm. In both cases, the focal spots are placed 25 mm from each other, which corresponds to the relative positions of (-25 mm), (0 mm) and (25 mm) from the optical axis in case of the symmetrical structure and (-50 mm), (-25 mm) and (0 mm) for the asymmetrical one. The following naming convention is adopted hereafter: R – right focal spot (25 mm); M – middle, on-axis focal spot (0 mm); L – left focal spot (-25 mm); 2L – far left focal spot (-50 mm). The structures resemble a combination of the Fresnel lenses (circles and rings) and a diffraction grating (vertical lines). This specific, combined design is also transmitted to the focal plane. In both cases, not only the three focal spots can be observed, but also the interferometric pattern in the form of vertical lines. This effect is undesirable, but the maximal values of the interference pattern are in the order of a few per cent of focal spot intensities and therefore are insignificant. Some differences in the intensity of the focal spots can be seen as well. In the symmetrical case, the ratios of intensities of the side focal spots to the intensity of the central one are equal to 99% and 98%. The stronger non-uniformity is observed for the asymmetrical case, where these ratios are equal to 93% and 78% for the L and 2L focal spots, respectively. These values have been obtained after several iterations of the optimization algorithm. The presented DOEs are phase diffractive structures coded as a kinoform. Therefore, their diffractive efficiencies are equal to 100% and for the plane wave illumination the whole radiation illuminating the structure is redirected into the three designed focal spots (here, not taking into account attenuation of the structure and the reflection losses). Additionally, a numerical method based on suppressing the values of the higher peaks has been implemented in order to reduce the differences between the intensities of the focal spots.

## 3. Manufacturing of diffractive optical elements

Both designed structures have been fabricated using additive manufacturing from polyamide PA 12 [16]. The dependence of the absorption coefficient and refractive index of various types of 3D-printable PA12 material on the frequency of radiation is presented in Fig. 2. The absorption coefficient of the material used for manufacturing measured at 150 GHz is equal to 0.5 cm$^{-1}$. Selective laser sintering method [4,5] has been used which allows to reproduce details as small as 0.05 mm. However, it should be also underlined that very thin and high element parts cannot reach such resolution. They must have a thickness of around 0.8 mm to form a high enough structure detail. In order to obtain a file that can be later 3D printed, the calculated phase distribution must be transformed into a height map $h(x,y)$, according to Eq. (1):

The height map is formed by multiplying the phase value $\phi (x,y)$ by the design wavelength $\lambda$ and dividing it by 2$\pi$ and the difference of the refractive indices of the used material $n=1.59$ and the surrounding material, which is air. Thus, the maximal height of the structure corresponds to $2\pi$ phase delay and is equal to 3.39 mm. The substrate layer of 1 mm was added, which is necessary to provide sufficient stiffness of the structure. Therefore, the total maximal height of the structures was 4.39 mm. It corresponds to the average absorption losses of 9%. Furthermore, the Fresnel reflection losses for the structures are equal to 5%. Since both structures were designed to have 200 mm of diameter, they were manufactured in the same size. Photographs of the 3D-printed lenses can be seen in Fig. 3.

## 4. Experimental verification

Manufactured structures have been verified in the experimental setup presented in Fig. 4. It consists of the frequency multiplier (VDI WR-5.1) for radiation emission at 150 GHz and the Schottky diode for detection. Two-dimensional scans have been obtained by moving the detector and gathering the intensity at multiple points in the focal plane. The measurement pitch was equal to 1.5 mm in the wide scans (panels (a) in Fig. 5 and Fig. 6) and 1.0 mm in the single focal point scans.

Structures-under-test have been placed at a distance $z_e =$ 1320 mm from the emitter treated as quasi-point source. According to the thin lens equation $\frac {1}{z_e}+\frac {1}{z_d}=\frac {1}{f}$ it resulted in focusing at a distance of $z_d =$ 108 mm. The actual longitudinal placement of the detector has been tuned to obtain the maximal signal across all measured focal points and has shown to be equal to $z_{sym} =$ 111 mm and $z_{asym} =$ 106 mm for symmetric and asymmetric designs, respectively. The obtained results are presented in Fig. 5 and Fig. 6.

In the symmetrically designed structure (Fig. 5) the L, M and R focal spots have been detected at positions equal to (-28 mm), 0 mm and 28 mm, respectively (two of them out of the main optical axis). The 3 mm difference from the designed positions (-25 mm) and 25 mm is connected with the image scaling, proportional to the illumination with a divergent beam (not parallel) resulting in different focusing distance, thus different shifts along X axis as well. The focal spot maxima on both sides are equal to slightly more than half of the central spot intensity. These values are significantly lower than expected, which is most likely caused by the directivity of the detector used. It is symbolically shown in Fig. 4. A conical horn is attached to the detector, which allows to detect incoming radiation with the full 3 dB bandwidth equal to $14^{\circ }$. It can be easily calculated that in the described geometry the light illuminates the detector at the angle of about $14^{\circ }$ from the center of the structure and up to about $51^{\circ }$ from the far edge of the structure. It means that the light collected by the detector at each measurement point comes from a certain area of the lens resulting from a finite angle of view of the detector. When raster scanning is performed, it can be understood as the detector integrating the radiation shaped by the structure. However, at any point the detector does not collect the radiation coming from the whole structure. This fact combined with the Gaussian-shaped illumination explains the lower than expected values measured at the side focal spots. Nevertheless, all three focal spots have been detected in the desired positions with reasonable intensities and the concept of the multi-focal-spot diffractive lens for the THz spectral range has been verified.

Also in the asymmetrical case [illustrated in Fig. 6(a)] all focal points have been registered. Their positions are equal to M – 0 mm, L – (-26 mm) and 2L – (-54 mm), which once again comes from the scaling factor connected with the setup geometry. Maximal intensities of the 2L and L focal spots correspond to 0.40 and 0.44 of the maximal intensity of the central one, respectively. As predicted, the intensities of the off-axis focal spots are smaller, the differences, however, are more significant than in the theoretically simulated patterns. Here, also the directivity of the detector is crucial – even more for the 2L focal spot, where the incidence angle can be as high as about $56^{\circ }$. Its intensity is therefore suppressed even stronger. This issue can be partially solved by rotation of the directive detector in a way that the center of the lens lies on the axis of the acceptance cone. The results of detector rotation of an angle of about $27^{\circ }$ are displayed in Figs. 6(e), 6(f), 6(g) and 6(h). At this angle, the far-left focal spot is observed from the center of the lens. The improvement is twofold – firstly, the radiation coming from the central part of the lens is detected and secondly, the radiation from the bigger part of the lens is gathered. It is connected with the fact that the acceptance cone is now placed at the rotation angle of the detector and is tilted in relation to the lens plane. Thus, it collects the radiation from the elliptical shaped area of the lens. The relative intensities of the focal spots obtained in all three measurements are summarized in Table 1. They are also complemented with additional measurements performed with the Schottky diode with the silicon hemispherical lens as a detector. Focal spots theoretically favored by the angle of the detector have been underlined in each row. Several conclusions can be formulated on this basis. First of all, the strongest focal spot is always placed on the optical axis (M). Secondly, the intensity values of both focal spots moved 25 mm from the optical axis (L and R) in relation to the central one are higher in the symmetrical design, when measured with the horn antenna. These proportions are comparable for the first and fifth rows in Table 1 thanks to the preserved geometry of the verification setup (detector registering angle). In the symmetrical design, they are equal to 0.55 and 0.51, which is on average $20\%$ higher in comparison to 0.44 obtained in the asymmetrical case. This behavior has also been predicted in the numerical simulations, however the differences in intensity were significantly smaller. Moreover, it can be noticed that for the rotated detector the results are much more favorable for the shifted focal spots reaching relative intensities of 0.54 and 0.78 for the 2L and L shifts, respectively. Measurements of the maximal intensities of the particular focal spots have also been performed with different detector – the Schottky diode equipped with the hemispherical lens. Measurements were conducted for both structures for various detector angles. A hemispherical lens is less directive than the horn antenna, however the influence of the direction of the incoming radiation is also clearly observable. The smaller angular sensitivity of the hemispherical lens resulted in the higher relative intensity values of the side focal spots in all measurements. For the asymmetrical structure one can observe similar behavior of intensity ratios as in the case of measurements conducted with the horn antenna. The rotation of the detector only partially decreases dominance of the M focal spot. Interestingly enough, the highest relative intensity of the L focal spot do not appear for the $14^{\circ }$ angle, but for $27^{\circ }$ (though the maximal measured absolute intensity indeed occurred for this angle). It is most likely connected with the strongest diminishing of the M focal spot intensity during the measurements at an angle of $27^{\circ }$, which translates to the highest L/M ratio of the intensities. As one can expect the highest relative intensity of the 2L focal spot has been observed for the $27^{\circ }$ configuration. Additionally, a procedure with the tilted hemispherical lens detector has been applied for the symmetrical structure, including both incidence angles, corresponding to L and R focal spots. As shown in the second row, for the detector angle of $0^{\circ }$, the higher relative intensity distribution – equal to 0.69 – has been registered for the R focal spot. It indicates that this focal spot is stronger comparing to the L focal spot with intensity equal to 0.59. The same behavior has also been observed for the measurements with the horn antenna, though it was less significant (the relative intensities of 0.55 and 0.51, respectively). In the $-14^{\circ }$ angle configuration, favoring the R focal spot, the highest intensity still was registered for M focal spot, located on the main optical axis. Nevertheless, the relative intensity for the R spot reaches 0.71, while for the L focal spot it drops to 0.45. On the other hand, the measurements favoring L spot, in the $14^{\circ }$ angle configuration, have not shown any differences between the relative intensities of L and M focal spots. The intensity of the R focal spot, however, strongly declined to 0.39. The comparison of the results obtained for different detector types and various angle configurations is extremely hard. One can observe some general tendencies, but any conclusions should be drawn very carefully. It shows beyond a reasonable doubt that the measurement technique and even the detector used are crucial for the reliability of the results. Some investigation in this topic has already been performed, proving the initial theses. Still, this is definitely a direction for the future work and there is plenty of space for the further improvement in the field of characterization of off-axis optical elements for the THz spectrum.

Nevertheless, the domination of the central (on-axis) focal spot in the experimental evaluation is observed in all detection geometries. It is therefore crucial to weaken the intensity of this focal spot in relation to all shifted ones. Below we propose a solution based on neural network training algorithms to equalize the focal spot intensities using a novel optimization method.

## 5. Optimization with neural networks

An algorithm for optimization of diffractive structures can be also realized with neural networks [17–21]. The iterative nature of net training shows visible similarities to the iterative algorithm described above (based on G-S algorithm). When the structure of the net realizes the propagation of light at a certain distance, the performed calculations can be understood as a forward propagation of the optical field. On the other hand, the error back-propagation mechanism used in the training process resembles the backward propagation of light in G-S algorithm. In our approach, the light propagation is also simulated with the convolution method. However, the implementation is realized by the neural network framework. An additional advantage is that the mathematical operation of convolution is commonly applied in neural networks for image processing and its usage is therefore well described. Ultimately, the numerical method of the light propagation remains the same, but the optimization method changes. The NN-based propagation and the optimization algorithm has been implemented in the Wolfram Language. The in-depth description of the algorithms, the source code and the adequate examples can be found in the Wolfram Community blog post. Structures analogous to the ones shown in the previous part of the article have also been simulated with the described neural network framework. The comparison of the results obtained with not-optimized and optimized structures using neural network convolution method are shown in Fig. 7.

First of all, it should be noted that the previously described optimization algorithm allowed to reduce differences between the values of intensity in particular focal spots. Moreover, the iterative algorithm slightly mitigated the interference pattern, manifested by the parallel lines in the focal plane. However, this pattern is still visible (even in the experimental results in the top images in Fig. 5 and Fig. 6). In both cases, all three focal spots are also noticeably elongated along the mentioned lines (vertical axis). This behavior results from the shape of the structure. In all simulations, the structures have been illuminated with the circular Gaussian beam. The circular aperture transforms into a Bessel-like point in the focal plane. The calculated intensity distribution has also clearly visible vertical lines laid over the sum of three diffractive Fresnel lenses – a single propagation method in Fig. 7(a). Therefore, the shape of the focal spots is defined by a cross-section of the Bessel function with a vertical line. The iterative G-S based optimization process provided small changes in the intensity details, but the overall shape of the structure (the phase delay map) and the resulting intensity distribution have been not modified significantly. However, the results of the neural network-based optimizations are much more interesting. The first striking issue is related to the shape of the phase delay map, which is totally different from the ones based on classical numerical simulations. The NN-based phase delay map does not resemble any of the diffractive shapes that could be predicted by human intuition. Nevertheless, the desired focusing pattern is obtained, as presented in the upper right (c) panel in Fig. 7. It can be easily seen that the background around the focal spots is much smaller and uniform. The focal spots are also more circular (no horizontal elongation), but still their maximal intensity values are slightly differing (they are not equal). Though, we believe that it could be solved with a proper parametrization of the training process and such goal is included in our future research scope.

## 6. Conclusions

Three-focal-spot diffractive lenses in the symmetrical and most importantly in the asymmetrical designs have been successfully developed and verified. Highly non-symmetric distribution of the focal spots requires special design and measuring techniques. Application of the iterative optimization algorithms allowed to diminish the unwanted focal spots entirely. However, the issues of the non-equal intensity of the spots and the non-uniform background have not been solved completely. As justified above, the focal spot placed on the optical axis usually tends to have a larger intensity value (tends to be stronger). Therefore, we propose a machine learning enhanced optimization algorithm based on the neural network design, which allows to almost completely get rid of the non-uniformity of the background as well as partially reduce the differences between the intensities of the focal spots. Moreover, neural networks could be utilized to improve the second generation of structures, which is one of our future goals. If the measured result would differ from the simulated one (for any reason), an additional layer in the neural network could be implemented to emulate this deviation. As a result, DOEs of the next generation should be able to shape the radiation more closely to the designed distribution. Another challenge to overcome is the measuring process. The directivity of the detector plays an important role in the case of structures with high numerical aperture or having large angles of deflection of light. Different detectors and measuring techniques are tested to obtain the results with the highest possible compliance with reality, which is still not satisfactory. The main drawback is the field of view of the detector, which is mostly too small to gather the whole off-axis radiation redirected by the optical structure.

## Funding

Narodowe Centrum Badań i Rozwoju (LIDER /11/0036 / L-9/17 / NCBR / 2018).

## Acknowledgments

The Authors would like to thank Institute of Optoelectronics at Military University of Technology for enabling experimental verification of samples. The Authors would like to thank Orteh Company for providing LS 6.0 software used for the design and modelling of diffractive optical elements, which is accessible in the Laboratory of Optical Information Processing at the Faculty of Physics in Warsaw University of Technology.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **A. Siemion, “Terahertz diffractive optics–smart control over radiation,” J. Infrared, Millimeter, Terahertz Waves **40**(5), 477–499 (2019). [CrossRef]

**2. **D. C. O’Shea, T. J. Suleski, A. D. Kathman, and D. W. Prather, * Diffractive optics: design, fabrication, and test*, vol. 62 (SPIE press, 2004).

**3. **E. Castro-Camus, M. Koch, and A. I. Hernandez-Serrano, “Additive manufacture of photonic components for the terahertz band,” J. Appl. Phys. **127**(21), 210901 (2020). [CrossRef]

**4. **U. M. Dilberoglu, B. Gharehpapagh, U. Yaman, and M. Dolen, “The role of additive manufacturing in the era of industry 4.0,” Procedia Manuf. **11**, 545–554 (2017). [CrossRef]

**5. **Z.-X. Low, Y. T. Chua, B. M. Ray, D. Mattia, I. S. Metcalfe, and D. A. Patterson, “Perspective on 3d printing of separation membranes and comparison to related unconventional fabrication techniques,” J. Membr. Sci. **523**, 596–613 (2017). [CrossRef]

**6. **K. Liebert, M. Rachon, A. Siemion, J. Suszek, D. But, W. Knap, and M. Sypek, “Thz beam shaper realizing fan-out patterns,” J. Infrared, Millimeter, Terahertz Waves **38**(8), 1019–1030 (2017). [CrossRef]

**7. **H. Dammann and K. Görtler, “High-efficiency in-line multiple imaging by means of multiple phase holograms,” Opt. Commun. **3**(5), 312–315 (1971). [CrossRef]

**8. **U. Krackhardt and N. Streibl, “Design of dammann-gratings for array generation,” Opt. Commun. **74**(1-2), 31–36 (1989). [CrossRef]

**9. **V. X. Yang, N. Munce, J. Pekar, M. L. Gordon, S. Lo, N. E. Marcon, B. C. Wilson, and I. A. Vitkin, “Micromachined array tip for multifocus fiber-based optical coherence tomography,” Opt. Lett. **29**(15), 1754–1756 (2004). [CrossRef]

**10. **J. B. Perraud, A. F. Obaton, J. Bou-Sleiman, B. Recur, H. Balacey, F. Darracq, J.-P. Guillet, and P. Mounaix, “Terahertz imaging and tomography as efficient instruments for testing polymer additive manufacturing objects,” Appl. Opt. **55**(13), 3462–3467 (2016). [CrossRef]

**11. **M. Sypek, “Light propagation in the fresnel region. new numerical approach,” Opt. Commun. **116**(1-3), 43–48 (1995). [CrossRef]

**12. **R. W. Gerchberg, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik **35**, 237–246 (1972).

**13. **R. G. Dorsch, A. W. Lohmann, and S. Sinzinger, “Fresnel ping-pong algorithm for two-plane computer-generated hologram display,” Appl. Opt. **33**(5), 869–875 (1994). [CrossRef]

**14. **A. Siemion, P. Komorowski, M. Surma, I. Ducin, P. Sobotka, M. Walczakowski, and E. Czerwińska, “Terahertz diffractive structures for compact in-reflection inspection setup,” Opt. Express **28**(1), 715–723 (2020). [CrossRef]

**15. **T. D. Ngo, A. Kashani, G. Imbalzano, K. T. Nguyen, and D. Hui, “Additive manufacturing (3d printing): A review of materials, methods, applications and challenges,” Composites, Part B **143**, 172–196 (2018). [CrossRef]

**16. **S. Busch, M. Weidenbach, M. Fey, F. Schäfer, T. Probst, and M. Koch, “Optical properties of 3d printable plastics in the thz regime and their application for 3d printed thz optics,” J. Infrared, Millimeter, Terahertz Waves **35**(12), 993–997 (2014). [CrossRef]

**17. **Y. Luo, D. Mengu, N. T. Yardimci, Y. Rivenson, M. Veli, M. Jarrahi, and A. Ozcan, “Design of task-specific optical systems using broadband diffractive neural networks,” Light: Sci. Appl. **8**(1), 112 (2019). [CrossRef]

**18. **X. Lin, Y. Rivenson, N. T. Yardimci, M. Veli, Y. Luo, M. Jarrahi, and A. Ozcan, “All-optical machine learning using diffractive deep neural networks,” Science **361**(6406), 1004–1008 (2018). [CrossRef]

**19. **L. Lu, L. Zhu, Q. Zhang, B. Zhu, Q. Yao, M. Yu, H. Niu, M. Dong, G. Zhong, and Z. Zeng, “Miniaturized diffraction grating design and processing for deep neural network,” IEEE Photonics Technol. Lett. **31**(24), 1952–1955 (2019). [CrossRef]

**20. **T. Yan, J. Wu, T. Zhou, H. Xie, F. Xu, J. Fan, L. Fang, X. Lin, and Q. Dai, “Fourier-space diffractive deep neural network,” Phys. Rev. Lett. **123**(2), 023901 (2019). [CrossRef]

**21. **C. Qian, X. Lin, X. Lin, J. Xu, Y. Sun, E. Li, B. Zhang, and H. Chen, “Performing optical logic operations by a diffractive neural network,” Light: Sci. Appl. **9**(1), 59 (2020). [CrossRef]