Subwavelength Diffractive Optical Elements for Generation of Terahertz Coherent Beams with Pre-Given Polarization State

Coherent terahertz beams with radial polarization of the 1st, 2nd, and 3rd orders have been generated with the use of silicon subwavelength diffractive optical elements (DOEs). Silicon elements were fabricated by a technology similar to the technology used before for the fabrication of DOEs forming laser terahertz beams with pre-given mode content. The beam of the terahertz Novosibirsk Free Electron Laser was used as the illuminating beam. The experimental results are in good agreement with the results of the computer simulation.

Structuring radiation from new sources in the terahertz range, including high-power ones, such as free-electron lasers (FELs) [41], requires optical elements designed with consideration of the features of such radiation, including wavelength and high-power density [42][43][44][45]. A good overview of the recently invented terahertz optical structures based on diffraction design is presented in [46]. The fabricated diffractive optical elements (DOEs) were used to focus [42,43,[47][48][49] and split [45,50,51] the terahertz laser beam, as well as to control the transverse-mode composition of the beam [52][53][54]. In particular, silicon binary elements were used to transform the illuminating beam of a high-power freeelectron terahertz laser into the Hermite-Gaussian, Laguerre-Gaussian, and Bessel singlemode beams [52,54]. However, such DOEs were used just to change the transverse mode composition without alteration of the polarization state of the illuminating beam. Note that some relevant applications of laser radiation require simultaneous control of the transversemode composition and polarization state of the beam [3,4,6,13,20,21]. There are well-known works on the polarization transformation of radiation in the terahertz range based on metaldielectric metasurfaces [55][56][57]. However, all-dielectric metasurfaces [58][59][60] are preferred because they are chemically inert and are not subject to oxidation.
Previously [59,60], the authors designed, fabricated, and examined a meta-axicon (axicon with a subwavelength period) for converting linearly polarized terahertz radiation into a second-order radially polarized beam. In this paper, we present new simulation and experimental results for the conversion of linearly polarized terahertz radiation into n TE e f f = Qn 2 1 + (1 − Q)n 2 2 1/2 (1) n TM e f f = Qn −2 where n 1 is the refractive index of the first medium, n 2 is the refractive index of the second medium, n TE e f f is the ordinary refractive index corresponding to the direction parallel to the layers of the structure, n TM e f f is the extraordinary refractive index corresponding to the direction perpendicular to the layers of the structure. Q = d 1 d 1 +d 2 is a fill factor (d 1 is the thickness of the silicon grating ridge and d 2 is the distance between adjacent grating ridges).
As practice shows, Formula (3) is not sufficiently accurate. Full-vector numerical calculation is required to determine the optimal height of the relief.
An important advantage of optical elements based on subwavelength gratings is the ability to change the directions of the slow and fast axes of the crystal by controlling the orientation of the grooves of the subwavelength grating. Thus, it is possible to create polarizing optical elements that convert the input linearly polarized radiation into beams with cylindrical polarization of various orders.
Let us consider the process of transmission through the subwavelength grating in Jones notation. Jones vector of the falling light has a view: Jones matrix of the half-wave element has a view: M λ/2 = cos 2ϕ sin 2ϕ sin 2ϕ − cos 2ϕ (5) where ϕ is an angle between the half-wave plate axis and the x-axis.
In the simplest case, when the second component of the Jones vector is equal to 0, the matrix (5) works like a rotation matrix that rotates incident polarization at an angle 2ϕ. Moreover, a subwavelength grating can have curved grooves causing different angles of rotation at different points of the grating. In that way, we can modulate a falling linear polarization and create a beam with spatially modulated polarization. Radially polarized cylindrical vector beam with a topological order of n has the following form: here θ is a polar angle.
To create the beam (6), we can modify an incident linearly polarized beam with a subwavelength grating that has a matrix: Comparing (7) and (5), we can conclude that the angle ϕ between the x-axis and fast axis of the grating has the following form: In general, the binary height of an element is determined by the formula: where (r, ϕ) are the polar coordinates, sign() is the sign function and f (r, ϕ) is the phase of the grating. The phase function of the element forming the beam (6) for m = 1 and 3 has the form [64]: For m = 2, the phase function of the element will have the form: In (10) and (11) d 0 is a constant that determines the period of the grating, ϕ 0 which is an angle between the beam orientation direction and horizontal axis x. Figure 1 shows a general view of fast and slow axes for the generation of differentorder radially polarized cylindrical vector beams. polarization and create a beam with spatially modulated polarization.
Radially polarized cylindrical vector beam with a topological order of n has the fo lowing form: here θ is a polar angle.
To create the beam (6), we can modify an incident linearly polarized beam with subwavelength grating that has a matrix: (7) and (5), we can conclude that the angle ϕ between the x-axis an fast axis of the grating has the following form: In (10) and (11) 0 d is a constant that determines the period of the grating, 0 ϕ whic is an angle between the beam orientation direction and horizontal axis x. Figure 1 shows a general view of fast and slow axes for the generation of differen order radially polarized cylindrical vector beams. Thus, we can design subwavelength gratings with curved grooves to generate cylindrical vector beams of various orders.

Design, Simulation, and Fabrication of Subwavelength Diffractive Optical Elements
The subwavelength optical element should have a height that corresponds to a halfwave plate. However, a subwavelength grating is not completely equivalent to a half-wave plate. Therefore, we use numerical simulation to find the height of the subwavelength grating that provides the best quality of the formed beam.
To find the optimal etching depth of a subwavelength grating, we consider an element of the order n = 2. As can be seen from Figure 1, this element has the form of an axicon and represents equidistant concentric annular ridges of a subwavelength grating.
The scheme of calculation is shown in Figure 2. The calculated area has the form of a cylinder. The subwavelength grating is located at the lower part of the domain; linearly polarized light passes through the subwavelength grating from bottom to up. Electric field amplitude distributions at different cross-sections are also shown in Figure 2. Thus, we can design subwavelength gratings with curved grooves to generate cylin drical vector beams of various orders.

Design, Simulation, and Fabrication of Subwavelength Diffractive Optical Elements
The subwavelength optical element should have a height that corresponds to a hal wave plate. However, a subwavelength grating is not completely equivalent to a hal wave plate. Therefore, we use numerical simulation to find the height of the subwave length grating that provides the best quality of the formed beam.
To find the optimal etching depth of a subwavelength grating, we consider an ele ment of the order n = 2. As can be seen from Figure 1, this element has the form of a axicon and represents equidistant concentric annular ridges of a subwavelength grating The scheme of calculation is shown in Figure 2. The calculated area has the form of cylinder. The subwavelength grating is located at the lower part of the domain; linearl polarized light passes through the subwavelength grating from bottom to up. Electric fiel amplitude distributions at different cross-sections are also shown in Figure 2. Let us consider the distribution of the electric field amplitude in the cross-section o the computational domain at different heights of the subwavelength gratings (Table 1 The first column shows the full amplitude of the electric field, the second column contain the x-component, and the third column shows the y-component. Let us consider the distribution of the electric field amplitude in the cross-section of the computational domain at different heights of the subwavelength gratings ( Table 1). The first column shows the full amplitude of the electric field, the second column contains the x-component, and the third column shows the y-component.
We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other.  We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is  We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is  We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is 50 We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is  We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. Figure 3 shows the dependence of maxima values of the x-component of the horizontal and vertical pairs, as well as the values of the maximum of the y-component on the height of the relief of the subwavelength grating.
The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is We formulate a criterion for the quality of the beam in the following way. Table 1 shows that horizontal and vertical polarizations form distribution patterns with four local maxima. Moreover, the distribution of the x-component has horizontally and vertically arranged pairs of maxima, and the distribution of the y-component has maxima located diagonally. Moreover, the values of the horizontally located maxima of the x-component differ from the values of the maxima located vertically. And the values of all four maxima of the distribution of the y-components are equal to each other. The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is The amplitudes shown in Figure 3 have absolute values and are not normalized. According to the dynamics of these graphs, it can be seen that initially, most of the energy is contained in the X-component. This is expected because the original field is X-polarized. Increasing the height of the relief to 37 microns allows us to slightly increase the energy in the Y-component by reducing it in the X-component. However, this leaves an asymmetry in the structure of the X-components (see the first row of Table 1). Therefore, we considered a further increase in the relief height to 50 microns, where, firstly, the intersection of three graphs is observed, and, secondly, a symmetrical structure in both transverse components (see the second row of Table 1). We believe that this situation corresponds to the formation of a second-order radial polarization. Unfortunately, in this case, a significant part of the energy is lost, which is scattered on a diffraction structure with high relief. The subwavelength elements have been designed by methods based on the rigorou light theory [64]. The following DOE parameters were chosen: the aperture diameter D 50 mm, discretization step s = 10 μm, and wavelength λ = 141 μm. Figure 4a-c shows th calculated binary subwavelength microrelief of meta-axicons for generating terahert beams with radial polarization of the first, second, and third orders, respectively. Also the meta-axicons add a focusing phase to the beam (NA = 0.3). In neighboring ring-shape Fresnel zones, subwavelength grating ridges are perpendicular to each other that provid a focusing phase in the output beam. Figure 4d-f presents the pre-given transverse distr bution of the beams (red color for horizontal polarization and green color for vertical po larization). The model under consideration has the following form. Linearly polarized light fall on the element, as shown in Figure 2. The meta-axicon forms a radially polarized beam o  Figure 3 shows that the three lines intersect at a relief height of 50 microns. We will choose this height for the following manufacturing of the element.
The subwavelength elements have been designed by methods based on the rigorous light theory [64]. The following DOE parameters were chosen: the aperture diameter D = 50 mm, discretization step s = 10 µm, and wavelength λ = 141 µm. Figure 4a-c shows the calculated binary subwavelength microrelief of meta-axicons for generating terahertz beams with radial polarization of the first, second, and third orders, respectively. Also, the meta-axicons add a focusing phase to the beam (NA = 0.3). In neighboring ringshaped Fresnel zones, subwavelength grating ridges are perpendicular to each other that provide a focusing phase in the output beam. Figure 4d-f presents the pre-given transverse distribution of the beams (red color for horizontal polarization and green color for vertical polarization).  The subwavelength elements have been designed by methods based on the rigor light theory [64]. The following DOE parameters were chosen: the aperture diameter 50 mm, discretization step s = 10 μm, and wavelength λ = 141 μm. Figure 4a-c shows calculated binary subwavelength microrelief of meta-axicons for generating terahe beams with radial polarization of the first, second, and third orders, respectively. A the meta-axicons add a focusing phase to the beam (NA = 0.3). In neighboring ring-sha Fresnel zones, subwavelength grating ridges are perpendicular to each other that prov a focusing phase in the output beam. Figure 4d-f presents the pre-given transverse dis bution of the beams (red color for horizontal polarization and green color for vertical larization).  Table 2 presents the results of the computer simulation of field distributions im diately behind the meta-axicons MAx1 (Figure 4a), MAx2 (Figure 4b), and MAx3 (Fig  4c) with added polarization analyzer rotated to the appropriate angle. The radius of domain was 1.5 mm.
The model under consideration has the following form. Linearly polarized light f  Table 2 presents the results of the computer simulation of field distributions immediately behind the meta-axicons MAx1 (Figure 4a), MAx2 (Figure 4b), and MAx3 (Figure 4c) with added polarization analyzer rotated to the appropriate angle. The size of the domain was 1.5 mm. Table 2. Results of computer simulation of field distributions immediately behind the meta-axicons that form radially polarized beams with the order of 1, 2, and 3. Beams pass through an analyzer rotating clockwise; the incident beam is polarized along the x-axis.

Analyzer Rotation Angle (Deg)
Radial Polarization Order n = 1 n = 2 n = 3 0 sector structure in accordance with the rotation angle of the polarization analyzer. Field distribution has a multi-ring structure because of the axicon-type focusing structure used in the metasurface. Table 2. Results of computer simulation of field distributions immediately behind the meta-axicons that form radially polarized beams with the order of 1, 2, and 3. Beams pass through an analyzer rotating clockwise; the incident beam is polarized along the x-axis. sector structure in accordance with the rotation angle of the polarization analyzer. Field distribution has a multi-ring structure because of the axicon-type focusing structure used in the metasurface. Table 2. Results of computer simulation of field distributions immediately behind the meta-axicons that form radially polarized beams with the order of 1, 2, and 3. Beams pass through an analyzer rotating clockwise; the incident beam is polarized along the x-axis. sector structure in accordance with the rotation angle of the polarization analyzer. Field distribution has a multi-ring structure because of the axicon-type focusing structure used in the metasurface. Table 2. Results of computer simulation of field distributions immediately behind the meta-axicons that form radially polarized beams with the order of 1, 2, and 3. Beams pass through an analyzer rotating clockwise; the incident beam is polarized along the x-axis.  Results of computer simulation of field distributions immediately behind the meta-axicons that form radially polarized beams with the order of 1, 2, and 3. Beams pass through an analyzer rotating clockwise; the incident beam is polarized along the x-axis.  Results of computer simulation of field distributions immediately behind the meta-axicons that form radially polarized beams with the order of 1, 2, and 3. Beams pass through an analyzer rotating clockwise; the incident beam is polarized along the x-axis.  Results of computer simulation of field distributions immediately behind the meta-axicons that form radially polarized beams with the order of 1, 2, and 3. Beams pass through an analyzer rotating clockwise; the incident beam is polarized along the x-axis.    Figure 5 shows the results of simulating the focusing of the formed beam in the focal plane of the meta-axicons. The amplitude distributions have a special structure consisting of 2n spots, which is consistent with the results of the article [62]. In scalar theory, a ring should be in focus, but in the article [62], it was shown that the presence of a longitudinal field component significantly affects the focal distribution, which distorts the picture of full intensity. The designed subwavelength elements (Figure 4a-c) were realized using the lithography technology used in [59,60] for the fabrication of terahertz subwavelength axicon. Previously, a similar technology based on the Bosch process [65] was used for the fabrication of diffractive optical elements to form laser beams with a pre-given orbital angular  Figure 5 shows the results of simulating the focusing of the formed beam in the focal plane of the meta-axicons. The amplitude distributions have a special structure consisting of 2n spots, which is consistent with the results of the article [62]. In scalar theory, a ring should be in focus, but in the article [62], it was shown that the presence of a longitudinal field component significantly affects the focal distribution, which distorts the picture of full intensity. The designed subwavelength elements (Figure 4a-c) were realized using the lithography technology used in [59,60] for the fabrication of terahertz subwavelength axicon. Previously, a similar technology based on the Bosch process [65] was used for the fabrication of diffractive optical elements to form laser beams with a pre-given orbital angular  Figure 5 shows the results of simulating the focusing of the formed beam in the focal plane of the meta-axicons. The amplitude distributions have a special structure consisting of 2n spots, which is consistent with the results of the article [62]. In scalar theory, a ring should be in focus, but in the article [62], it was shown that the presence of a longitudinal field component significantly affects the focal distribution, which distorts the picture of full intensity. The designed subwavelength elements (Figure 4a-c) were realized using the lithography technology used in [59,60] for the fabrication of terahertz subwavelength axicon. Previously, a similar technology based on the Bosch process [65] was used for the fabrication of diffractive optical elements to form laser beams with a pre-given orbital angular

Analyzer Rotation
The model under consideration has the following form. Linearly polarized light falls on the element, as shown in Figure 2. The meta-axicon forms a radially polarized beam of order n with a focusing phase NA = 0.3. Next, the beam passes through an analyzing polarizer. The orientation angle of the polarizer axis relative to the horizontal axis is given in the left column. After passing the polarization analyzer, maxima are allocated for 2n sectors in accordance with the order of polarization n. Table 2 shows the rotation of the sector structure in accordance with the rotation angle of the polarization analyzer. Field distribution has a multi-ring structure because of the axicon-type focusing structure used in the metasurface. Figure 5 shows the results of simulating the focusing of the formed beam in the focal plane of the meta-axicons. The amplitude distributions have a special structure consisting of 2n spots, which is consistent with the results of the article [62]. In scalar theory, a ring should be in focus, but in the article [62], it was shown that the presence of a longitudinal field component significantly affects the focal distribution, which distorts the picture of full intensity. Figure 5 shows the results of simulating the focusing of the formed beam in the focal plane of the meta-axicons. The amplitude distributions have a special structure consisting of 2n spots, which is consistent with the results of the article [62]. In scalar theory, a ring should be in focus, but in the article [62], it was shown that the presence of a longitudinal field component significantly affects the focal distribution, which distorts the picture of full intensity. The designed subwavelength elements (Figure 4a-c) were realized using the lithography technology used in [59,60] for the fabrication of terahertz subwavelength axicon. Previously, a similar technology based on the Bosch process [65] was used for the fabrication of diffractive optical elements to form laser beams with a pre-given orbital angular moment [54]. An SEM image of a realized element microrelief (Figure 4b) is shown in Figure 6. The designed subwavelength elements (Figure 4a-c) were realized using the lithography technology used in [59,60] for the fabrication of terahertz subwavelength axicon. Previously, a similar technology based on the Bosch process [65] was used for the fabrication of diffractive optical elements to form laser beams with a pre-given orbital angular moment [54]. An SEM image of a realized element microrelief (Figure 4a) is shown in Figure 6. Note that the simulation results (Table 2 and Figure 5) showed an incomplete polarization transformation, as well as the complexity of the analysis from patterns in the near diffraction zone. Therefore, to filter parasitic components, an optical scheme with a Fourier analyzer was used in the experimental study.
The following causes of incomplete polarization transformation using the meta-axicon can be distinguished. Firstly, these are Fresnel reflections due to the high refractive index (n1 = 3.42) of the meta-axicon material. Secondly, there is a non-uniform reflection due to the dependence of the refractive index on the orientation of the grating grooves. Thirdly, the thickness of the subwavelength grating also significantly affects the transmission since the grating works as a thin film. Fourth, when designing the element, an ideal averaging of the refractive indices of the lattice material and air was assumed, but the Note that the simulation results (Table 2 and Figure 5) showed an incomplete polarization transformation, as well as the complexity of the analysis from patterns in the near diffraction zone. Therefore, to filter parasitic components, an optical scheme with a Fourier analyzer was used in the experimental study.
The following causes of incomplete polarization transformation using the meta-axicon can be distinguished. Firstly, these are Fresnel reflections due to the high refractive index (n 1 = 3.42) of the meta-axicon material. Secondly, there is a non-uniform reflection due to the dependence of the refractive index on the orientation of the grating grooves. Thirdly, the thickness of the subwavelength grating also significantly affects the transmission since the grating works as a thin film. Fourth, when designing the element, an ideal averaging of the refractive indices of the lattice material and air was assumed, but the experiment was carried out in a more rigorous model, which showed a lower efficiency of the element than the ideal model.

Experimental Investigation of Subwavelength Diffractive Optical Elements
The fabricated subwavelength meta-axicons (an image of one of the axicons taken with an electron microscope is shown in Figure 6) were investigated employing the terahertz radiation of the Novosibirsk free electron laser at Budker Institute of Nuclear Physics of the Siberian Branch of the Russian Academy of Sciences [41]. Laser radiation emerging as a continuous stream of 100-ps pulses with a repetition rate of 5.6 MHz was tuned to a wavelength of 141 µm. The spectral width of the radiation was about 1 µm. The Gaussian beam, x-polarized with a wire polarizer, incident on a meta-axicon, as shown in Figure 7. Cross-section of the beam passed through the meta-axicon was recorded ( Figure 7a) by a pyroelectric camera Pyrocam IV with a matrix of 320 × 320 pixels (the size of one element was 80 µm). The total image size was 25.6 × 25.6 mm 2 . Since the camera is not sensitive to radiation polarization, in this configuration, we observed the distribution of the total beam intensity, regardless of its local polarization. For the exploration of the local polarization of the beam, an analyzer (Figure 7b) was introduced into the system. Rotating the analyzer, we were able to study the distribution of local polarization in the beam.

PEER REVIEW
10 of 16 the total beam intensity, regardless of its local polarization. For the exploration of the local polarization of the beam, an analyzer (Figure 7b) was introduced into the system. Rotating the analyzer, we were able to study the distribution of local polarization in the beam.
In both configurations (with and without an analyzer), a polypropylene kinoform lens 80 mm in diameter with a focal length of 75 mm could be added to the optical system, as shown in Figure 7c. In this case, the recording part of the optical system turns into a Fourier analyzer, which makes it possible to study the beam spectrum behind the metaaxicon in the space of transverse wave numbers. The large diameter of the lens makes it possible to register both positive and negative orders of laser beam diffraction by the meta-axicon. A photograph of the experimental setup is shown in Figure 7d. In the case of a diffractive axicon, the expected intensity distribution in the focal plane is a ring [66]. The experiments were carried out using the axicons of the first and third orders. The purpose of the experiments was to compare the experimentally measured distributions of the local polarization of the beam after passing through the axicon illuminated by linearly polarized radiation and the polarization distribution obtained by numerical simulations. In both configurations (with and without an analyzer), a polypropylene kinoform lens 80 mm in diameter with a focal length of 75 mm could be added to the optical system, as shown in Figure 7c. In this case, the recording part of the optical system turns into a Fourier analyzer, which makes it possible to study the beam spectrum behind the meta-axicon in the space of transverse wave numbers. The large diameter of the lens makes it possible to register both positive and negative orders of laser beam diffraction by the meta-axicon. A photograph of the experimental setup is shown in Figure 7d. In the case of a diffractive axicon, the expected intensity distribution in the focal plane is a ring [66].
The experiments were carried out using the axicons of the first and third orders. The purpose of the experiments was to compare the experimentally measured distributions of the local polarization of the beam after passing through the axicon illuminated by linearly polarized radiation and the polarization distribution obtained by numerical simulations. First experiments were carried out using the optical scheme of Figure 7b. They showed that behind both studied axicons, the Gaussian beam transforms into beams with a singularity near the optical axis, as predicted by the calculations (see Figure 4). The introduction of the analyzer into the optical system ( Figure 7c) made it possible to see that for axicons of the first and third order, the local direction of polarization corresponds to the theoretically expected one.
For the MAx1 axicon at the beam periphery, as expected according to Figure 4a, the observed intensity distribution depended on the azimuthal angle α as cos 2 α. In the case of the MAx3 axicon, the intensity depended on the angle as cos 2 (3α). A selection of frames recorded with Pyrocam IV for this axicon at several positions of the analyzer is shown in Figure 8. The scheme of rotation of the intensity maxima for this case is easy to understand using Figure 9. and Table 3. When the analyzer is rotated by α, the observed six annual sectors are rotated by β = α/3 (see Table 3).
Sensors 2023, 23, x FOR PEER REVIEW 11 of 16 understand using Figure 9. and Table 3. When the analyzer is rotated by α , the observed six annual sectors are rotated by Table 3).   Installing a lens in the optical scheme allowed us to obtain more detailed information Figure 8. Experimentally observed intensity distribution in the beam that passed through diffractive element MAx3 and the analyzer with polarization direction α; β-observed image rotation angle shown with yellow arrows (compare with Figure 9).  Installing a lens in the optical scheme allowed us to obtain more detailed information about the formed beams. For both axicons, in the absence of an analyzer, a uniform ring is observed in the focal plane of the lens (Figure 10a,d). This means that the beams formed are, in essence, a kind of Bessel beams of first and third orders, and they are a superposition of conically converging plane waves produced by the axicons. Since most of the beam energy is concentrated in the ring, we may assume that the diffraction efficiency is rather high. Quantitative measurements of diffraction efficiency will be carried out later. If an analyzer is installed between the element and the chamber, then one can observe changes in the intensity distribution associated with the polarization properties of the beam. In particular, when installing the analyzer, the intensity of the ring in the focal plane will be modulated in azimuth, and 2n segments will be observed, where n is the order of the axicon. In the case of the MAx1 axicon, two half-rings are observed (Figure 11b). When the analyzer is rotated by α = 40 0 , the pattern rotates by the same angle (Figure 11c).
For the MAx3 axicon, when the analyzer is rotated by an angle α , the observed six annual sectors are rotated by an angle β . For the analyzer set in the y-direction, we observed the segments in 30 ±  comparing to the vertical (Figure 10e). This case corresponds to Figure 9 ( 90 α =  ). Rotation of the analyzer by 90  leads to the pattern shown in Figure 10f, which is schematically illustrated in Figure 9 ( 0 α =  ). In support of the above, at the end of the section, we present one more set of images of the beam crosssection in the focal plane of the lens, obtained by rotating the analyzer for the MAx3 axicon ( Figure 11). If an analyzer is installed between the element and the chamber, then one can observe changes in the intensity distribution associated with the polarization properties of the beam. In particular, when installing the analyzer, the intensity of the ring in the focal plane will be modulated in azimuth, and 2n segments will be observed, where n is the order of the axicon. In the case of the MAx1 axicon, two half-rings are observed (Figure 10b). When the analyzer is rotated by α = 40 • , the pattern rotates by the same angle (Figure 10c).
For the MAx3 axicon, when the analyzer is rotated by an angle α, the observed six annual sectors are rotated by an angle β. For the analyzer set in the y-direction, we observed the segments in ±30 • comparing to the vertical (Figure 10e). This case corresponds to Figure 9 (|α| = 90 • ). Rotation of the analyzer by 90 • leads to the pattern shown in Figure 10f, which is schematically illustrated in Figure 9 (α = 0 • ). In support of the above, at the end of the section, we present one more set of images of the beam cross-section in the focal plane of the lens, obtained by rotating the analyzer for the MAx3 axicon ( Figure 11). observed the segments in 30 ±  comparing to the vertical (Figure 10e). This cas corresponds to Figure 9 ( 90 α =  ). Rotation of the analyzer by 90  leads to the patter shown in Figure 10f, which is schematically illustrated in Figure 9 ( 0 α =  ). In support o the above, at the end of the section, we present one more set of images of the beam cross section in the focal plane of the lens, obtained by rotating the analyzer for the MAx3 axico ( Figure 11).

Conclusions
Silicon subwavelength diffractive optical elements for the generation of terahertz coherent beams with radial polarization of the first and third orders have been designed and fabricated by reactive ion etching. Methods of rigorous light theory have been used for the design of subwavelength microrelief. The fabricated elements were investigated by methods of computer simulation and natural experiments. The beam of the terahertz Novosibirsk Free Electron Laser was used as the illuminating beam (wavelength λ = 141 µm). The experimental results are in good agreement with the results of the computer simulation. It was experimentally and numerically shown that the used approach allows the generation of terahertz coherent beams with the possibility to control transverse mode content and polarization state simultaneously. This possibility is crucial for such applications as lidars [3,4], terahertz telecommunication systems with MDM [3], remote control systems for fundamental research, and so forth. In [67], it was discussed that a relatively high ratio between terahertz range wavelength and optical material structuring resolution opens an opportunity to create optical elements forming terahertz coherent fields with pre-given intensity, phase, mode content as well as polarization state.

Data Availability Statement:
The data that support the findings of this study are available from the authors upon reasonable request.