Neural Network-Based Design of Two-Focal-Spot Terahertz Diffractive Optical Elements

This article presents a convolutional neural network approach for the design and optimization of single-input-multiple-output (SIMO) structures for the sub-terahertz spectral range (140 GHz). Two SIMO structures with two output channels have been designed using the proposed neural network approach and an iterative algorithm as a reference. Both structures have been manufactured by means of fused deposition modeling 3-D printing technique and verified experimentally. A new method of 3-D modeling of the designed phase maps has been developed and applied to manufacture unintuitive structures optimized with neural networks.

Neural Network-Based Design of Two-Focal-Spot Terahertz Diffractive Optical Elements Paweł Komorowski , Marta Kurowska , Mateusz Kaluza , Patrycja Czerwińska , Przemysław Zagrajek , and Agnieszka Siemion Abstract-This article presents a convolutional neural network approach for the design and optimization of single-input-multipleoutput (SIMO) structures for the sub-terahertz spectral range (140 GHz).Two SIMO structures with two output channels have been designed using the proposed neural network approach and an iterative algorithm as a reference.Both structures have been manufactured by means of fused deposition modeling 3-D printing technique and verified experimentally.A new method of 3-D modeling of the designed phase maps has been developed and applied to manufacture unintuitive structures optimized with neural networks.

I. INTRODUCTION
I N THE last few decades, terahertz (THz) radiation has al- ready found numerous applications in various areas, such as telecommunication, nondestructive testing, detection of dangerous objects, or medical diagnostics [1], [2], [3], [4], [5], [6], [7], [8].Despite the remarkable progress already made, many crucial technologies, including the generation, detection, and manipulation of the THz waves, still require further development.This work focuses on the latter aspect, namely the methods of forming and redirecting the THz beams into the given shapes and patterns.One of the crucial points of the application of THz waves is the next generation of wireless telecommunication.Future telecommunication systems will require the distribution of signals to end-user devices with constantly increasing amounts of transferred data [9].Wi-Fi technology is currently commonly used for short-range wireless data transfer, reaching its limits in the data capacity of links.The solution is switching to higher frequencies in the THz radiation band [10].Moreover, techniques to increase the data capacity of telecommunication links, known from other systems, should also be applied to the THz-based transmission band.The crucial technology used in modern wireless data transmission systems is known as multiple-input-multiple-output (MIMO) [11].A variant of this method based on a combination of a multiple-input-singleoutput (MISO) setup, a common optical data transmission link, and a single-input-multiple-output (SIMO) setup could be proposed for THz communication.Such a system would operate in free space and guarantee multiplexing of radiation in terms of time, frequency, spatial signal location, or polarization [12].
Structures similar to SIMO have already been considered for lower frequencies, especially microwave K and Q bands.Diffractive elements, known also as reflectarrays [13], [14], [15] or transmitarrays [16], [17], [18] have been proposed for various applications including space technology and telecommunication.Although dealing with substantially different challenges due to different wavelengths of radiation considered, such arrays also share a few similarities with the diffractive optical elements (DOEs) analyzed in this article.In both cases, the spatial distribution of the phase retardation (phase delay map) is introduced to the wavefront of radiation to alter its shape.Interest in the structures realizing multibeam patterns in much more matured microwave technology indicates the need for the development of such elements also for the emerging THz technology.
DOEs can be considered as an interesting solution to form THz beams into arbitrary optical field patterns [19].Such elements can be designed using different, even very sophisticated methods.For the simplest patterns, such as a single focal point or focal line, the appropriate analytical equations can be used to describe the shape of a phase mask.More complex radiation shapes require, however, more advanced methods based on the simulated propagation between structure and focal plane in an iterative manner.A novel idea is to use a neural network (NN) for the emulation of radiation propagation and use the NN training process for the optimization of the DOE.
It is worth noting that, for the sub-THz and THz frequencies, DOEs can often be manufactured with relatively simple and cost-efficient 3-D printing methods, while much more advanced techniques, such as photolithography or electron-beam lithography, are necessary for infrared or visible radiation ranges [20].Additive manufacturing techniques such as fused deposition modeling (FDM), selective laser sintering (SLS), and digital light processing (DLP) can be used in the DOE fabrication process for the THz radiation range [21].These methods use polymer materials in different forms in the fabrication process.Recent research indicates that some materials used in the FDM method, such as polypropylene (PP), high impact polystyrene (HIPS), or cyclic olefin copolymer (COC), have a very low absorption coefficient in the THz range [22], [23].Thus, they can be used during the manufacturing process of phase passive optical components required in the case of a THz SIMO element.Moreover, 3-D printing techniques provide high-quality and fabrication resolution of the DOEs in the THz radiation range.Therefore, 3-D-printed diffractive elements are an excellent solution for fast and inexpensive prototyping of THz beamforming systems.To adjust this technique to fast-varying phase maps obtained through the NN algorithm, a new 3-D modeling method has been developed and applied to manufacture these structures.
In this work, we present SIMO structures that allow spatial multiplexing of incoming THz radiation, designed using an iterative algorithm and a neural network approach.The proposed NN-based optimization method was used for the first time in the research described in this article and two previous works [24], [25].It allowed for the design of unique structures with nonintuitive phase distributions, unattainable by other methods.Using NNs to design or optimize diffraction structures is a new concept.Nevertheless, in recent years, the use of NNs has gained increasing interest [26], [27], [28].Diffractive deep neural networks (D 2 NNs) have so far been used to classify objects and numbers [29], [30], [31], perform image segmentation [32] and logical operations [33], or shape and steer radiation beams [34], [35], [36], [37], [38], [39].An unquestionable advantage and distinguishing feature of D 2 NNs is the ability to perform calculations at the speed of light.After setting up the optical system corresponding to the designed network, the calculations are carried out by displaying the intensity distribution corresponding to the considered image (e.g., the object to be categorized) on the system input.The light propagates along the D 2 NN and diffracts on its successive elements, finally focusing it in the appropriate position in the output plane.This location depends on the initially displayed image and varies for different input images.One should name two fundamental differences between the network described in this article and the solutions discussed above.First, it is the number of hidden layers, in this case, equal to 1.The proposed algorithm aims to design single-plane diffractive optical structures for THz spectral range.Another significant change is the use of the convolutional method of light propagation instead of fully connected layers, based directly on the Huygens principle.Convolutional neural networks (CNNs), in turn, are commonly used in machine vision and object recognition [40], which makes these types of networks a relatively well known and refined method.
The idea standing behind the network presented in this article is to combine unique optical properties of D 2 NN with computational efficiency and refined framework of CNNs.The proposed NN uses a convolution layer, but in a different way than typical CNNs [41].In the CNNs, usually the convolution kernels undergo changes in the training process, which allows for optimizing, for example, prediction models.A completely different situation occurs in the proposed algorithm, where the diffraction integral defines the kernel.In this approach, an additional layer describing the phase map of the optimized structure is added to the network and the weights of this layer are changed in the training process.Nevertheless, the mathematical operation of convolution itself does not change, and from this perspective, the proposed network is computationally similar to CNNs.Therefore, it is a hybrid network that uses a convolution operation and realizes the diffraction of light.To the authors' knowledge, the convolutional diffractive neural network (CDNN) has not been used or described in the literature before.
In this article, the CDNN approach is compared with a betterknown Gerchberg-Saxton (GS) iterative algorithm in the task of design and optimization of a two-focal-spot diffractive lens.As it will be shown, this novel approach allowed to obtain completely new phase maps describing new DOEs, which realize intensity patterns very similar to the DOEs optimized with the established GS method.This fact alone is already interesting from the perspective of the theory of diffraction, however, the NN method offers much more possibilities.The most interesting prospect is a diffractive element optimized for multiple frequencies (which in general is impossible for typical diffractive optics).However, the same design method, based on neural network training on the set with varying parameters, can also be applied, for example, for varying illuminating intensity patterns, the direction of the incoming radiation, its divergence, or even nonflat surfaces of the optical elements.

A. Design
Computer-generated holography is a useful and smart tool when it comes to designing DOEs, especially those designed to work in the THz regime.It allows one to form the radiation into almost arbitrary shapes.Choosing the best design method that will suit the desired performance of the structure is crucial.Designing DOEs with multiple focal points requires more precise and sublimed techniques than designing a classic Fresnel lens, especially when the lens is meant to work in the off-axis regime.Only the simplest DOEs with a single on-axis focal point can be designed using theoretical equation (depending on the configuration -in on-axis or off-axis approaches) or repropagation while designing a structure with well-shaped shifted off-axis focal point is already challenging [42].The case of multifocal-spot lenses is practically infeasible with these classic methods.The most common solution here is the usage of iterative procedures, based on the GS algorithm [43].Even such a solution can be not enough to obtain the desired distributions of the focal spots with uniform intensity and circular shape, which brings the necessity of introducing modifications that count for additional aberrations.Another novel approach is to use an NN for the DOE design and optimization process.The proposed method utilizes an NN to simulate radiation propagation and supervised training for the design and optimization of the structure.These two design methods are described and compared in this section.
The iterative GS algorithm simulates the propagation of radiation between the hologram and image planes and overwrites the amplitude of the optical field in those planes with the corresponding amplitude distributions without changing its phase.Fig. 1.Scheme of the GS algorithm.In the central rectangle, the complex optical field is propagated N times between the input and output planes.The input (Gaussian) and target (cross) amplitudes are forced in the input and output planes, respectively.After optimization, the phase distribution in the input plane defines the optimized element, while the square of the amplitude distribution in the output plane shows the predicted intensity distribution behind the element.
The scheme of such an algorithm in the form of a flowchart is presented in Fig. 1.The complex optical field is visualized as a color image, where the brightness corresponds to the amplitude while the change of color from blue through green and yellow to red indicates the phase shift (from 0 to 2π).In the hologram plane, the amplitude of the incident beam, in this case, the Gaussian distribution, is forced.In the image plane, the forced amplitude corresponds to the desired amplitude distribution (a cross in Fig. 1 has been chosen to emphasize the influence of the optimization algorithm on the output; in the actual simulations, the target image consisted of two focal spots).In both planes, the phase distributions are calculated and do not change according to predefined distributions.The relation existing between these two planes is realized by forward and backward propagation.The propagation is based on a convolution theorem [44], [45], where the complex optical field is convolved with the impulse response function, as in where U 0 and U 1 denote complex optical fields in the input and output planes, respectively, while (x,y) are the Cartesian coordinates in both these planes.The h function is known as the impulse response function of the free space and is expressed according to where k = 2π λ is the wavenumber, z is the propagation distance, λ is the wavelength, and i is the imaginary unit.
The optimization starts in the input plane, where the input amplitude is fed to the algorithm without any information about the phase distribution (the diffractive structure).This optical field is then propagated forward to the output plane.Here, the absolute value of the complex optical field (real amplitude) is overwritten with the amplitude distribution of the target image (a cross in Fig. 1).From this point, the complex amplitude describing the optical field is propagated backward to the original input plane, where similarly the amplitude of the optical field is overwritten with the input amplitude.In both steps of this algorithm, the phase distribution is not altered.This allows coding information about the target amplitude distribution in the phase distribution in the input plane.The whole process is repeated several times, which allows to optimize the phase map of the designed element.The gain dependence on the number of iterations for the more complex multiplane modification of the GS algorithm was described by Makowski [46], [47].The optimal number of iterations of the algorithm slightly depends on the complexity of the forced amplitude.However, in the cases investigated in this article, only a few iterations are necessary to obtain the optimized phase delay map of the designed element.
From a theoretical and descriptive point of view, the radiation propagation method used in the NN algorithm does not differ from the one used in the GS algorithm.In both cases, the square matrix of pixels containing information about the complex amplitude of radiation is discretely convolved with another matrix of the same size and shape containing sampled impulse response functions for the free space [the h function as in (2)].The main difference between these two approaches lies in the optimization algorithm.In the NN scheme, the information about the difference between the obtained and target output amplitudes is propagated backward, while in the GS scheme, only the target amplitude is backpropagated.This more complex information paves the way for more advanced optimization algorithms [48], based on gradient descent or adaptive moment estimation.On the other hand, such algorithms require additional layers of calculations for reshaping the input and output, determining the value of the loss function, and defining the structure with trainable parameters.The flowchart of the layers in the NN scheme is shown in Fig. 2.
The NN can be described with subsequent layers joined with different links between particular nodes.The input layer (red nodes) defines the distribution of the optical field illuminating the structure (U 0 ) as well as propagation parameters, such as wavelength or distance.These parameters are used to describe the convolution core [the impulse response function h, defined by (2)].After the input layer, there are hidden layers gathered into three groups.The first hidden layer (green node) is the only trainable layer of the network.It consists of the matrix of phase Fig. 2. Scheme of the convolutional diffractive neural network (CDNN).The nodes symbolically describe subsequent layers of the network.Arrows show the flow of the data within the network.The input data in the form of the amplitude of illuminating radiation and propagation parameters are given to the network (red nodes).In the green node the input field is multiplied with the phase map defining the diffractive structure (these are trainable weights of the network).Next, four blue nodes perform four real convolutions between real and imaginary parts of the optical field (U) and impulse response function (h), which are then combined back into a complex amplitude in the dark blue node.Finally, the output layer (violet node) provides the output from the network in the form of the real amplitude of the optical field in the focal plane.retardation values defined for every pixel of the simulated optical structure.The mathematical operation performed by this layer is an element-by-element multiplication of the input field (U 0 ) with the phase retardation matrix described as a complex amplitude Values of the phase matrix (φ) are modified in the network training process to guarantee the best matching between the obtained and target output amplitudes.The next hidden layers denoted with blue nodes realize the crucial mathematical operation -convolution of the optical field with an impulse response function.The convolution of two complex matrices has been described using four real convolutions (4) which are then combined back into a single complex matrix that describes the optical field of the output layer (dark blue node) Therefore, in Fig. 2 four parallel convolution layers are shown (blue nodes) followed by a summation layer (light blue node), which adds and subtracts convolution results into a single complex matrix according to (4).
The square of the modulus of this complex amplitude describes the intensity pattern observed in the output plane.The last layer (violet node) provides the output of the algorithm in the form of a matrix of real values of field amplitude.This value is used later on to compare the calculated intensity distribution with the desired one and determine the value of the loss function.As described above, an additional layer, which defines the phase map of the optimized structure, is being changed in the training process.In this approach, the weights of the convolution kernels are constant and defined according to the real and imaginary parts of (1).Weights of the second layer (trainable layer in Fig. 2) are optimized in the training process and define the phase map of the designed structure.The back-propagation algorithm optimizes the structure by minimizing the loss function (mean squared error), defined as where i and j are horizontal and vertical indices of the calculation matrix of size n × n, U 1 and U t are the calculated and target amplitude distributions, respectively.The adaptive moment estimation (ADAM [49]) method has been used in this case.To guarantee that the algorithm stops at the right moment, the best-fit result is kept in the memory and overwritten with the better result when such is obtained.The optimization stops when the best result has not changed for a given number of iterations-usually a few tens or hundreds.The end condition has been chosen empirically and depends on the parameters of the network training process, such as the learning rate.In the training process, a single-element training set is used.The goal of optimization is to redirect the radiation to the given intensity distribution exactly.Typically, such an approach would result in overfitting and is generally avoided.However, in this approach, this is a desired outcome, as a single best-fit structure is being developed to realize a particular diffractive task.
The fundamental data format used in a proposed NN is a 3-D matrix of the size 2 × n × n.This size comes from the separate notation of the real and imaginary parts of the complex amplitude, described as a square matrix with n pixels on the edge.The dual matrix is separated into matrices describing only real or imaginary values before the convolution operations.The outputs of the convolution (n × n matrices) are then once again combined into a dual matrix (see Fig. 2).
It is worth noticing that the process of NN training naturally uses very similar concepts as iterative optimization algorithms known from computer-generated holography.In the diffractive NN, the flow of data from the input to the output corresponds to the propagation of the radiation at a given distance.On the other hand, from the algorithmic point of view, the back-propagation of error is analogous to the propagation at a negative distance.However, one should emphasize that the NN backpropagates more complete information-not only about the desired output but also about the difference between the desired and obtained results.

B. Modeling
Two SIMO elements with two output channels have been designed using GS and NN algorithms.The design wavelength (DWL) of both structures was equal to 2.14 mm which corresponds to the frequency of 140 GHz.The choice of DWL was dictated by the lower edge of the frequency band, generated by the available emitter.In this way, the smoothest manufacturing Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. is ensured, as the wavelength is the longest.It must be noted, however, that there are no contraindications to using this method also for higher frequencies.The focal spots of DOEs had been placed symmetrically with respect to the main optical axis of the system-25 mm to the right and 25 mm to the left of it.The focal length was equal to 700 mm.Such geometrical parameters ensured that the radiation was redirected at a small angle (±2 • ), minimizing the influence of the directivity of the detector used for the experimental verification of the structures [50].DOEs had been designed and simulated in square pixel matrices with pixel side lengths equal to 0.9 mm.This sampling was dictated by the nozzle size of the FDM 3-D printer used for manufacturing.In this way, two lines of material (0.45 mm each) were extruded into a single pixel of the structure, which turned out to be optimal during the testing of the additive manufacturing process.The GS lens had been trimmed to a circle with 150 mm in diameter, while the NN structure was shaped in a square with a side equal to 115 mm (128 × 128 matrix).These shape differences resulted from the implementation details of particular design algorithms.
The CDNN has been trained under different conditions with varying initial weights, learning rates, input amplitude distributions, and target images.The initial weights and shape of the input amplitude have shown minimal impact on the results and the calculation time.On the other hand, the influence of the learning rate and shape of the target intensity distributions is significant.The learning curves obtained for different training parameters are shown in Fig. 3.
The dependence of the loss curves on the learning rate (α) is presented in Fig. 3(a).The remaining parameters of the ADAM method have been set to β 1 = 0.9, β 2 = 0.999 and = 10 −5 .In the first few hundred rounds, the higher learning rate guarantees faster learning and better-fit result.However, after a few thousand rounds with high learning rates (α = 0.8 and α = 0.4), the algorithm falls into unstable oscillations and fails to converge.On the other hand, smaller learning rates result in slower learning at the beginning of the training process but allow for obtaining lower loss values.It must be noted, though, that the smallest learning rate does not necessarily return the best possible results.First of all, too low value of the learning rate unnecessarily stretches the calculation time, and secondly, for a fixed optimization time or a number of rounds, it will not be able to converge to the loss values obtained with a slightly higher learning rate.In the presented evaluation the learning rate of α = 0.1 performed best and therefore it has been used in further calculations.Fig. 3(b) presents the loss curves obtained with the fixed learning rate α = 0.1 and changing target distributions.Two focal spots in the output plane have been designed as Gaussian bells with 3σ width equal to the Airy diffraction limit.Subsequently, this width has been modified by a factor k to model sub-diffraction-limit spots (k < 1) and wider ones (k > 1), according to where σ is the standard deviation of the Gaussian distribution, k is the scaling factor, r A is the radius of an Airy disk, λ is the wavelength, f is the focal distance and d is the diameter of the aperture.
As can be seen in Fig. 3(b), for a fixed number of learning rounds the best results have been obtained for the focal spots slightly wider than the diffraction limit (k ∈ (1.25, 2.0)).It can be intuitively understood because the wider focal spots offer more nonzero points in the target image, which can be used by the algorithm for more precise adjustment.This effect has its limits, as for much wider spots (k = 4) the functioning of the DOE moves from focusing to collimation of the beam.Much more important, however, are the results for sub-diffraction-limit focal spots.It can be observed that the smaller the spot, the worse the optimization is.Ultimately, for k = 0.25 the algorithm gets stuck in an unoptimized position and does not improve regardless of the training time.It shows the importance of selecting proper, physically feasible target images for the optimization process.The unscaled focal spots (k = 1) have been chosen for optimization as a compromise between the efficiency of focusing and trainability.
The phase maps obtained with GS and NN methods are shown in Fig. 4 together with the corresponding simulated intensity distributions in the focal plane.In both cases, the target of the optimization has been defined as two diffraction-limited, Gaussian-shaped spots, separated by 50 mm.The spots have been placed symmetrically with respect to the main optical axis of the setup.The differences in the shape describing the phase retardation map of both structures can be easily noticed.The GS structure partially resembles a combination of two Fresnel lenses that have different parameters, joined in alternating stripes.Though it must be noted, that none of these theoretical lenses would focus the radiation outside of the optical axis.They are both symmetrical with respect to the center of the element and therefore, each of them separately should focus the radiation on the optical axis.However, their combination allows redirecting the radiation into off-axis focal spots, which can be seen in Fig. 4(c).The effects of the optimization algorithm can also be seen in the discontinuities between particular zones and stripes of the structure.These small adjustments can be clearly seen in the form of bright and dark pixels near the borders of the circular zones.Ideal Fresnel lenses defined with the theoretical equation would have sharp edges between particular zones (borders of white and black areas).
On the other hand, the structure obtained through the NN optimization design does not resemble any recognizable shape.From a human point of view, it is unintuitive and chaotic.The changes in the phase retardation introduced by this DOE are mostly discontinuous.Phase maps exhibit frequent 0-2π drops or alternating high and low values.Nevertheless, such an element reshapes the illuminating radiation in the same way as the previously described, more intuitive GS structure [Fig.4(c) vs. Fig.4(d)].In both cases, two focal spots with uniform intensity, placed at the desired positions and circularly shaped, have been obtained.In the NN design case, small cross-like disturbances can be observed from the aperture influence of the focal spots, which is connected with the shape of the aperture in the simulations, square versus circular, used in the GS case.The circular aperture transforms into the Bessel function in the Fourier space, while the square aperture results in the sinc function in the x and y dimensions.The aperture can be adjusted to the same shape; however, due to the computational memory limitation,

C. Manufacturing
The size of the details of the DOEs is related to the wavelength.For the DWL of 2.14 mm the hundreds of micrometer resolution is perfectly sufficient and can be met by relatively simple and widely available 3-D printing techniques, such as FDM, SLS, or DLP.These methods use different forms of polymer material: filament, powder, and resin, respectively.They allow for low-cost, fast prototyping, and manufacturing of complex structures [19], [51], which is crucial in the field of optics.Moreover, the THz time-domain spectroscopy shows that some polymer materials used in the FDM method have desired optical properties and can be used for manufacturing of phase kinoforms [52].Based on the results of our previous studies [19], [53], these materials are characterized by a very flat dispersion curve with a refractive index close to 1.5-1.6.At the same time, some of them do not introduce unwanted radiation attenuation due to the low absorption coefficient and loss tangent in the lower part of the THz frequencies spectrum.Fig. 5 illustrates the optical properties of COC material compared with two polymers commonly used in FDM 3-D printing as follows: 1) polyamide (PA12) and 2) polylactic acid (PLA).The data presented are based on measurements performed with THz time-domain spectroscopy.The COC has a significantly lower absorption coefficient and loss tangent than two other polymers.Low radiation damping for COC material occurs in the entire spectrum presented in Fig. 5.For DWL of 2.14 mm (140 GHz) absorption coefficient of COC is equal to 0.061 cm −1 and the refractive index approaches 1.51.Thus, the COC has been chosen to manufacture SIMO structures.It is worth noticing, that both the design and manufacturing methods are well suited for the wide part of the sub-THz and THz bands.The DWL corresponding to 140 GHz has been chosen due to high manufacturing accuracy in relation to the wavelength as well as the availability of an efficient source of radiation.However, the optical elements manufactured with 3-D printing methods can be applied even up to 1 THz and stay competitive with other design and manufacturing methods.The manufactured elements are shown in Fig. 6.The FDM 3-D printing technology guarantees high-quality manufacturing of phase DOEs with a vertical print resolution of 100 μm.The shape of each structure was determined on the basis of grayscale bitmaps, where the brightness represented by the pixel values corresponded to the phase change introduced by the structures.The thickness of the structure is defined as in

III. RESULTS AND DISCUSSION
The manufactured elements have been verified in the experimental setup, shown in Fig. 7.The Schottky diode-based frequency multiplier has been used as a source of THz radiation.It emits a linearly polarized, strongly monochromatic beam of tunable frequency.The baseband frequency after multiplication by a factor of 18 radiates in the WR-5.1 band.It means that the source emits electromagnetic waves in the 140-220 GHz range.For our experimental verification, the lowest possible frequency was chosen.The source generates at this conditions radiation of 0.5 mW power.The attached conical horn antenna propagates with 21 dBi of gain at 13 • of full 3 dB beamwidth.The radiation was then collimated into a quasi-plane wave with a parabolic mirror and redirected to illuminate the designed DOEs.A Schottky diode placed in the WR-5.1 waveguide equipped with a conical horn antenna has been used as a detector.The detector was mounted on the 3-D motorized stage, which allowed it to gather raster scans in arbitrary planes behind the tested DOE.In addition, a lock-in amplifier was used to improve the signal-to-noise ratio (SNR) of the observed intensity pattern.The SNRs obtained during the measurements reached 55 for the GS structure and 40 for the NN structure.It comes from the fact that the maximal intensity registered for the GS structure was higher than for the one designed using NN, while the noise level remained roughly on the same level.
The scans measured behind the GS structure in the planes perpendicular to the optical axis (xy plane in Fig. 7) and parallel to it (xz plane in Fig. 7   along the optical axis, visible in the left part of Fig. 8(b), is expected and related to the nature of diffraction when focal spots are formed.
Two analogous structures designed using the NN optimization algorithm and fabricated with interpolation and cuboid methods have also been experimentally verified.The results are presented in Fig. 9.
The differences in the registered intensity distributions for these two structures are tremendous.For the cuboid 3-Dmodeling method the desired two-focal-spot operation has been observed, while for the triangle 3-D modeling method the obtained intensity pattern has very few resemblances with the designed one.The xz scan admittedly reveals the longitudinal characteristics of two focal spots but they are hardly distinguishable from the background.Furthermore, the pattern observed in the xy plane consists of many randomly spaced spots.The improvement introduced by the cuboid 3-D-modeling technique is, therefore, undeniable, and the application of such a method is necessary for manufacturing structures designed with the NN.The DOE designed with the NN and manufactured with the cuboid method correctly redirected the THz radiation into two focal spots according to the project.This is a very important result, showing that the holograms realizing particular diffractive tasks do not necessarily have to be determined.There might be (and, as shown in this case, there are) alternative shapes of the phase maps, allowing the radiation to be redirected into very similar patterns.These shapes are often unintuitive and cannot be obtained with any other method.One can also suspect that such an approach could allow to reshape the radiation in a way that is unachievable by other methods.However, this hypothesis requires further research.

IV. CONCLUSION
The CDNN presented in this work allows the design and optimization of DOEs for the THz spectral range.Radiation propagation is simulated in the NN approach using a convolution method to calculate the diffraction integral and the network training process leads to the adjustment of the phase map of the structure to redirect the radiation as close to the defined intensity distribution as possible.
Two two-focal-spot SIMO structures have been designed using the proposed NN and an iterative algorithm as a reference.Both elements have been manufactured by means of FDM 3-D printing technology and verified experimentally.In both cases, the THz radiation has been redirected into two well-separated focal spots.In addition, a new method has been developed for the preparation of 3-D models from phase maps.It is based on the extrusion of cuboids instead of triangle interpolation and allows for the manufacturing of complex structures designed with NN-based algorithms.
The phase delay maps obtained through the NN optimization are unintuitive and do not resemble the patterns known from DOEs designed using other methods.Nevertheless, they allow to reshape the incoming radiation into intensity distributions very similar to the ones obtained with the classic structures.Algorithms typically used in optics and computer-generated holography are deterministic and always return very similar phase delay maps.What is new about the proposed NN-based algorithm is that it is able to return completely new phase delay maps and subsequently diffractive structures, realizing the same or very similar optical fields as the traditional algorithms.The NN algorithm goes beyond typical DOE design methods and offers new types of diffractive structures.It is especially compelling in the THz and sub-THz spectral range, where a huge need for new solutions for beam manipulation is observed.Moreover, manufacturing and proof testing of complex DOEs is relatively simple in this band, as the size of the structure details depends linearly on the wavelength.The diffractive structures proposed in this article are designed primarily for multiplexing in THz wireless data transmission systems.It must be noted, however, that the potential area of application of the DOEs designed using CDNN is much wider and includes any kind of solution, where THz beamforming is relevant.

Fig. 3 .
Fig. 3. Values of the loss function (mean squared error) obtained for NN training with different parameters.(a) Loss curves for initial learning rates ranging from α = 0.01 to α = 0.8.(b) Loss curves for the fixed learning rate (α = 0.01) and varying dimensions of the target focal spots.The values in the legend correspond to the linear scaling of the focal spots' diameters with k = 1.0 defined as a diffraction-limited spot [as in (6)].

Fig. 4 .
Fig. 4. Phase maps of the two-focal-spot lenses optimized with (a) GS and (b) neural network (NN) algorithms compared with the corresponding simulated intensity distributions in the focal plane (c) and (d), respectively.The black color in the phase maps denotes the relative 0 phase retardation (minimal height of the structure), while the white color correlates with the 2π phase retardation (maximal height).

Fig. 6 .
Fig. 6.Photographs of structures manufactured with the cuboid extrapolation technique and optimized using (a) Gerchberg-Saxton (GS) and (b) neural network (NN) algorithms.The cuboids have a base of 0.9 mm x 0.9 mm and heights ranging from 0 to 4 mm.

Fig. 7 .
Fig. 7. Experimental setup used for the verification of two-focal-spot lenses.The collimated THz beam illuminated the designed diffractive optical element (DOE), which redirected the radiation into two focal points.The detectors were rotated to point at the center of the structure.
) are shown in Fig. 8.The xy scan has been performed in the focal plane of the structure (adjusted manually for the highest intensity in the focal spots) and includes both off-axis focal spots as well as the empty space between them.The xz scan has been registered from the plane lying 125 mm before the focal plane to the plane lying 100 mm behind the focal plane.It incorporates two THz beams, forming designed focal spots in the middle and defocusing afterward.As can be seen, two distinct focal spots have been registered.Focusing of the THz radiation is clearly visible in the xy plane as well as in the xz plane [panels (a) and (b) in Fig. 8, respectively].The focal spots are placed 25 mm to the left and right of the optical axis.They are uniform in intensity and have the desired geometrical shape.In the xz plane, one can observe long focusing distances, which proves that the incoming radiation is indeed separated into two independent beams.The residual energy

Fig. 8 .
Fig. 8. Experimental results measured for the Gerchberg-Saxton (GS) structure showing intensity cross sections in the (a) xy plane and (b) xz plane.

Fig. 9 .
Fig. 9. Experimental results measured for the neural network (NN) structures manufactured with the cuboid method (top row) and the interpolation method (middle and bottom rows).(a) xy plane for the cuboid diffractive optical element (DOE).(b) xz plane for the cuboid DOE.(c) xy plane for the triangle DOE.(d) xz plane for the triangle DOE.(e) Wider area scan of xy plane for the triangle DOE (white rectangle shows the original scan).(f) Wider area scan of xz plane for the triangle DOE (green rectangle shows the original scan).