# Frequency Division Multiplexing of Terahertz Waves Realized by Diffractive Optical Elements
## 2 Materials and Methods
### Numerical Simulations  
The performance of the structures has been simulated using a convolution based light propagation algorithm [38] implemented in Wolfram Mathematica. It emulates the propagation of radiation by convolving the complex amplitude of the optical field with the impulse response function of the free space. For the two structures with shorter focal lengths (\(f=100\) mm and \(f=200\) mm), calculations have been performed on 4096 \(\times\) 4096 matrices of 0.117 um square pixels. For the other two structures, the matrix size has been increased to 8192 \(\times\) 8192 due to the higher perpendicular separation of the focal spots for longer propagation distances. Firstly, all structures have been verified for the design frequency of 300 GHz. Next, the frequency of the incoming radiation has been blue- and red-shifted with the step equal to 30 GHz. Simulated positions and dimensions of the focal spots for the two structures with shorter focal lengths are illustrated in the form of intensity distributions in Figure 2.  
These distributions present fragments of the focal planes cut in such a way that the optical axis is placed in the middle of the left edge. The 0 th order of diffraction can be seen in the form of dots and rings around the optical axis. For the DWL, it is nonexistent. However, the further the frequency of illuminating radiation is shifted, the stronger non-diffracted signal is observed. As expected, the longer wavelengths are focused farther from the optical axis into more distorted spots. Two effects contribute to this distortion. Firstly, the diffraction limit of light (defined as the diameter of Airy spot) depends linearly on the wavelength, hence the widening of the spots for longer wavelengths. Secondly, for the off-axis focusing, the circularity of the focal spot decreases approximately linearly with the tangent of the beam deflection angle [39], hence the horizontal elongation of the focal spots. It can be noticed that, in the case of the \(f=100\) mm structure, the spots obtained for the frequencies differing by 30 GHz are hardly separable. For the \(f=200\) mm structure,the separation is better, especially in the shorter wavelength division. Focal spots for longer wavelengths are still overlapping but can be easily distinguished. It shows that longer focal lengths should provide better separation of the signals of different frequencies, which has been verified for the \(f=300\) mm and \(f=400\) mm structures in Figure 3.  
Figure 1: Phase distributions of four diffractive structures with different focal lengths: (**a**) \(f=100\) mm; (**b**) \(f=200\) mm; (**c**) \(f=300\) mm; and (**d**) \(f=400\) mm. White color corresponds to \(2\pi\) phase shift, while a black to no phase shift is introduced by the structure.  
Figure 2: Calculated positions of the focal spots for structures with the focal lengths \(f=100\) mm (**top**) and \(f=200\) mm (**bottom**). Intensity distributions are shown with a main optical axis in the middle of left edge of the illustrated area. Colors denote the frequency of the illuminating radiation as in the legend.  
It can be noticed that, at these focal distances, the signals spaced by 30 GHz are clearly spatially separated. Moreover, radiation of higher frequency can be demultiplexed by the designed DOEs (up to 450 GHz, which corresponds to a 150 GHz shift from the DWL).