“Perfect” Terahertz Vortex Beams Formed Using Diffractive Axicons and Prospects for Excitation of Vortex Surface Plasmon Polaritons

: Transformation of a Bessel beam by a lens results in the formation of a “perfect” vortex beam (PVB) in the focal plane of the lens. The PVB has a single-ring cross-section and carries an orbital angular momentum (OAM) equal to the OAM of the “parent” beam. PVBs have numerous applications based on the assumption of their ideal ring-type structure. For instance, we proposed using terahertz PVBs to excite vortex surface plasmon polaritons propagating along cylindrical conductors and the creation of plasmon multiplex communication lines in the future (Comput. Opt. 2019, 43, 992). Recently, we demonstrated the formation of PVBs in the terahertz range using a Bessel beam produced using a spiral binary silicon axicon (Phys. Rev. A 2017, 96, 023846). It was shown that, in that case, the PVB was not annular, but was split into nested spiral segments, which was obviously a consequence of the method of Bessel beam generation. The search for methods of producing perfect beams with characteristics approaching theoretically possible ones is a topical task. Since for the terahertz range, there are no devices like spatial modulators of light in the visible range, the main method for controlling the mode composition of beams is the use of diffractive optical elements. In this work, we investigated the characteristics of perfect beams, the parent beams being quasi-Bessel beams created by three types of diffractive phase axicons made of high-resistivity silicon: binary, kinoform, and “holographic”. The amplitude-phase distributions of the ﬁeld in real perfect beams were calculated numerically in the approximation of the scalar diffraction theory. An analytical expression was obtained for the case of the binary axicon. It was shown that a distribution closest to an ideal vortex was obtained using a holographic axicon. The resulting distributions were compared with experimental and theoretical distributions of the evanescent ﬁeld of a plasmon near the gold–zinc sulﬁde–air surface at different thicknesses of the dielectric layer, and recommendations for experiments were given.


Introduction
Photon beams with an orbital angular momentum (OAM), or vortex beams, are well known and have been used in optics for 25 years. The history of their research and fields of application were described in a recent review [1]. More recently, the application of OAM beams in the terahertz range began. A review of methods for generating vortex beams in this spectral range can be found in [2]. One of the applications of vortex beams in nanooptics is the generation of concentrically converging surface plasmon polaritons (SPPs) on flat conducting surfaces. To excite them, the vortex beam illuminates narrow circular or spiral slits that are cut, for example, in a metal-insulator-metal sandwich [3] to excite an SPP on the metal layer on the opposite side. At a certain spatial-phase distribution of the electromagnetic field of plasmons, their orbital angular momentum relative to the axis of the system will be nonzero. Surface plasmons of this kind have been obtained experimentally in the visible range in many works (see, for example, [4][5][6] and reviews [7,8]).
It is possible, however, to excite another type of vortex plasmons, which propagate along an axisymmetric conductor and simultaneously rotate around its axis. In this case, the Poynting vector of a plasmon moves helically over the conductor, and, therefore, the plasmon carries an orbital angular momentum. In Ref. [9], we proposed using a combination of vortex surface plasmon polaritons propagating along a metal wire to create a multiplex plasmon transmission line, which is similar to using vortex beams for multiplex communication in free space [10]. Since the propagation length of plasmons in the visible range is of the order of 10 µm [11], in transmission lines of macroscopic dimensions, it is necessary to use terahertz or mid-infrared radiation. In [12], experiments performed at a wavelength of 141 µm demonstrated for the first time the excitation of vortex plasmons on an axisymmetric conductor and their propagation over a distance of 100 mm. To excite plasmons on the cylinder, higher-order vortex beams illuminated its end face (the end-fire coupling technique [13]). In this excitation method, the spatial distribution of the complex amplitude of the electric field of the vortex beam and its overlap with the distribution of the vortex plasmon field play an important role. The vortex plasmon field distribution is determined both by the properties of the conductor material and the wavelength of the incident radiation, as well as by the topological charge of the exciting beam.
When creating a multiplex line, it is necessary to simultaneously launch several vortex plasmons with different topological charges on a cylindrical conductor. Therefore, the beams exciting these plasmons must have equal diameters, and all wave energy must be concentrated in one annular beam. Theoretically, so-called perfect beams possess such properties [14,15]. Such beams are generated due to focusing by a lens of Bessel beams of different orders with the same transverse wave numbers [16]. In practice, however, "perfect beams" are far from perfectness. The reason is that, in reality, there are no ideal Bessel beams with an infinite cross-section and infinite energy. The properties of real Bessel beams depend on the methods of their generation, which have even a greater effect on the characteristics of perfect beams, as will be shown below.
In the visible range, Bessel beams are usually generated via formation of phase transparencies using commercially available spatial light modulators. No radiation modulators are available so far in the terahertz range, and phase diffractive elements [17] made of plastics [18] or silicon [19] are mostly used. This work is devoted to analytical and numerical study of the characteristics of ideal beams obtained using diffractive phase axicons. In Section 2, we describe the characteristics of three types of axicons that can be used to generate Bessel beams. In Section 3, we explain the experimental setup in which ideal beams will be used, present an analytical calculation of the electric field of an ideal beam for the case of a binary phase axicon, and perform a comparison with numerical calculations. In Section 4, ideal beams created by different types of axicons are compared. Section 5 is devoted to a discussion of the results and their use in experiments on excitation of vortex plasmons.

Bessel Beams
The electrical field of Bessel beams (BBs) is described by the expression where J (κr) is the Bessel function of the first kind of the order , k = (k 2 z + κ 2 ) 1/2 , and ϕ is the azimuth angle. Hence, they are (see, for instance, [2]) a superposition of plane waves: and the Fourier amplitudes of which are described by the following expression: In optics, the quantum number is called the topological charge. From the above expressions, one can see that the ideal Bessel wave has infinite transverse size, infinite energy, and, like a plane wave, cannot exist in reality. A superposition of Bessel waves, however, can be used for an analytic description of really existing radiation modes-for example, Laguerre-Gaussian beams.
As seen from Equations (1)-(4), some properties of Bessel beams are of interest in terms of use in real-life optical systems. The maximum beam intensity is in the vicinity of the beam optical axis; the intensity in the beam cross-section does not change with distance; higher-order beams have an orbital angular momentum (OAM). Due to these properties, the beams are "non-divergent". The intensity distribution in the beam cross-section is capable of "self-recovery" after the beam's passage through obstacles and during the beam's propagation in turbulent media [20]. The beams can transfer angular momentum to micromechanical systems and particles [21,22].

Methods of Forming Quasi-Bessel Beams
Equations (2)-(4) show how beams with properties close to those of ideal Bessel beams can be formed in real life. Here, we will call them quasi-Bessel beams (QBBs). To this end, it is necessary to use conically converging beams. A zero-order Bessel beam, for example, is formed behind an axicon (conical lens) illuminated by a beam with a plane wavefront. If a spiral phase plate is added to the axicon [13], then a higher-order Bessel beam can be formed (Figure 1a). The phase transmission function of such a hybrid element, which can be called a spiral axicon (SA), is described by the expression where α is the cone-apex angle, and n is the refractive index of the material. Such axicons were fabricated using a 3D printer to generate Bessel beams at a radiation frequency of 0.3 THz in [18]. Since we are interested in radiation frequencies that are an order of magnitude higher and greater radiation powers, such volumetric plastic elements cannot be used because of the strong absorption of radiation and thermal destruction of plastic (see Figure 8 in [14]). In additions, the elements have the drawback of the ability to operate only at one given wavelength. This problem can be solved by forming quasi-Bessel beams via the use of phase diffraction gratings with a given phase shift distribution, which are also called diffractive axicons [17]. Figure 1b-d presents three types of diffractive axicons. The phase function of a binary axicon (BA) is described by a stepwise profile where and p = 2π/κ is the axicon period.   This problem can be solved by forming quasi-Bessel beams via the use of phase diffraction gratings with a given phase shift distribution, which are also called diffractive axicons [17].
A kinoform axicon (KA), as seen in Figure 1c, differs from the previous one in the fact that within each period, the phase incursion decreases linearly along the radius by the value 2π ( ) where the function fix( ) z converts a real number into an integer, discarding the fractional part. Essentially, this axicon is a diffractive analogue of a refractive spiral axicon (SA) (see Equation (5)). Another axicon, which we, following [23], will call the holographic axicon (HA), has a profile that is step-wise in radius, but continuous in azimuth (Figure 1d). The width of the zones of this axicon, unlike the previous ones, is not constant, but is determined by the position of the zeros of the Bessel function [24]. In our notation, this function can be written as follows: This expression is derived from (5), but the spiral plate is divided into zones along the radius. The zones are between the zeroes of a Bessel function of a given order, in which the point of reference for the azimuth angle ϕ, which is an argument of the function ( , ) r Φ ϕ , is shifted by / Δϕ = π  for odd zones using the Heaviside function A kinoform axicon (KA), as seen in Figure 1c, differs from the previous one in the fact that within each period, the phase incursion decreases linearly along the radius by the value 2π where the function fix(z) converts a real number into an integer, discarding the fractional part. Essentially, this axicon is a diffractive analogue of a refractive spiral axicon (SA) (see Equation (5)).
Another axicon, which we, following [23], will call the holographic axicon (HA), has a profile that is step-wise in radius, but continuous in azimuth ( Figure 1d). The width of the zones of this axicon, unlike the previous ones, is not constant, but is determined by the position of the zeros of the Bessel function [24]. In our notation, this function can be written as follows: This expression is derived from (5), but the spiral plate is divided into zones along the radius. The zones are between the zeroes of a Bessel function of a given order, in which the point of reference for the azimuth angle ϕ, which is an argument of the function Φ(r, ϕ), is shifted by ∆ϕ = π/| | for odd zones using the Heaviside function For the sign of the beam topological charge to change, the transparency must be rotated through 180 degrees about the axis lying in the plane of the diffractive element and passing through its center.
All the elements described above make it possible to form QBBs with the same intensity distribution in the cross-section, which coincides well with the Bessel function, and a phase function that corresponds to the exponent in Equation (1), but exists only on a certain limited length (after which they collapse): where D is the diameter of the element. Hereinafter, we will assume that the angles θ k are small. In this case, beam propagation can be described in the approximation of the scalar theory of diffraction. The so-calculated intensity distribution over the beam cross-section for the case of = 3 is shown in Figure 1e. Hence, we can conclude that diffractive axicons can be used in any optical systems in which it is required to form Bessel beams whose cross-section practically does not differ from that of the ideal BBs. Later on, however, we will see that a more detailed analysis reveals important differences in the characteristics of these beams.
Let us now discuss the technical side of the issue. In the visible range, phase transparencies can be generated, as discussed in the introduction, using commercially available spatial light modulators (SLMs) based on liquid crystals, whereas in the long-wave range, it is necessary to make optical elements from plastic or other materials transparent in this range. In the mid-infrared and near-terahertz ranges, materials absorb more strongly, and making helicoidal axicons with a continuous or step-wise profile requires the use of complicated technologies (multistep photolithography or laser ablation [25]). In addition, high-power free-electron lasers and wavelength-tunable sources of mid-infrared and terahertz radiation require elements with high radiation resistance and that are capable of operating with the radiation wavelength tuned within several tens of a percent. Hence, subject to the above requirements, binary axicons ( Figure 1b) made of high-resistance silicon are the most suitable (at least in the first experiments). This work will mostly address these axicons. Bessel beams with an OAM (vortex beams) in the terahertz range were first obtained using such axicons [26].

Transformation of Quasi-Bessel Beams into "Perfect Vortex Beams"
As follows from expression (1), the diameters of the rings of Bessel beams grow with the topological charge . For some applications, for example, when radiation is input into dielectric fibers, this property is a disadvantage. In the case of Bessel beams, the problem can be solved by transforming them using a lens. From (1)-(4), it follows that, regardless of the magnitude and sign of the topological charge, if the beams have the same transverse wave number κ, they form an annular beam of the radius in the focal plane of the lens, the beam topological charge being equal to the charge of the initial beam. Such beams are usually called "perfect vortex beams" [14]. The amplitude-phase distribution in the real "perfect vortex beams", however, depends on the method used to create the "parent" quasi-Bessel beam. In [27], the authors showed that the radius and width of the ring of ideal beams obtained by focusing of Bessel-Gaussian beams increased with the topological charge and demonstrated that in experiments at a wavelength of 623 nm.

Optical Scheme and Coordinate Sysem
For the reader to understand the meaning of the calculations performed in this section, as well as the characteristics of axicons and wavelengths selected for numerical examples, we present here a diagram of the experimental setup ( Figure 2) for which these calculations are carried out. A linearly polarized Gaussian beam I = I 0 exp(2r 2 /w 2 ) from the Novosibirsk free electron laser (NovoFEL) with a plane wavefront (wavelength of 141 µm and mode radius w = 11 mm) enters the user station through two radiation-resistant tungsten wire polarizers (P), the characteristics of which are described in [28]. The first polarizer is to regulate the beam power, and the second polarizer provides vertical polarization of the beam. Then, the beam passes through the beam polarization converter, which forms a vector beam [29] with the radial direction of polarization. In the experiments, the vector beam is formed using a Mach-Zehnder interferometer [30] or a segmented half-wave plate [31]. The laser radiation is an infinite train of 100-ps pulses with a repetition rate of 5.6 MHz. The average power of the NovoFEL Gaussian beam in the experiments is 50-100 W. Even after 50% attenuation of the beam because of the Fresnel reflection on the axicon, the beam energy was sufficient to fire a paper sheet, which was also partially transparent for terahertz radiation, in several seconds ( Figure 3). Preliminary experiments showed that after the vector beam passed through the BA and the lens, an annular beam was formed at the lens focus, which retained its radial polarization.
polarizer is to regulate the beam power, and the second polarizer provides vertical polarization of the beam. Then, the beam passes through the beam polarization converter, which forms a vector beam [29] with the radial direction of polarization. In the experiments, the vector beam is formed using a Mach-Zehnder interferometer [30] or a segmented half-wave plate [31]. The laser radiation is an infinite train of 100-ps pulses with a repetition rate of 5.6 MHz. The average power of the NovoFEL Gaussian beam in the experiments is 50-100 W. Even after 50% attenuation of the beam because of the Fresnel reflection on the axicon, the beam energy was sufficient to fire a paper sheet, which was also partially transparent for terahertz radiation, in several seconds ( Figure 3). Preliminary experiments showed that after the vector beam passed through the BA and the lens, an annular beam was formed at the lens focus, which retained its radial polarization.

Vortex Beams behind Binary Axicon
In [26], the authors used binary phase axicons similar to that shown in Figure 1b to form Bessel beams of the first and second orders in the terahertz range. The cross-section of one of these beams, which was recorded using a microbolometer array detector, is shown in Figure 3. In [26], it was shown analytically that the electric field of a wave at a distance z′ behind the axicon can be written in the paraxial approximation as tungsten wire polarizers (P), the characteristics of which are described in [28]. The first polarizer is to regulate the beam power, and the second polarizer provides vertical polarization of the beam. Then, the beam passes through the beam polarization converter, which forms a vector beam [29] with the radial direction of polarization. In the experiments, the vector beam is formed using a Mach-Zehnder interferometer [30] or a segmented half-wave plate [31]. The laser radiation is an infinite train of 100-ps pulses with a repetition rate of 5.6 MHz. The average power of the NovoFEL Gaussian beam in the experiments is 50-100 W. Even after 50% attenuation of the beam because of the Fresnel reflection on the axicon, the beam energy was sufficient to fire a paper sheet, which was also partially transparent for terahertz radiation, in several seconds (Figure 3). Preliminary experiments showed that after the vector beam passed through the BA and the lens, an annular beam was formed at the lens focus, which retained its radial polarization.

Vortex Beams behind Binary Axicon
In [26], the authors used binary phase axicons similar to that shown in Figure 1b to form Bessel beams of the first and second orders in the terahertz range. The cross-section of one of these beams, which was recorded using a microbolometer array detector, is shown in Figure 3. In [26], it was shown analytically that the electric field of a wave at a distance z′ behind the axicon can be written in the paraxial approximation as

Vortex Beams behind Binary Axicon
In [26], the authors used binary phase axicons similar to that shown in Figure 1b to form Bessel beams of the first and second orders in the terahertz range. The cross-section of one of these beams, which was recorded using a microbolometer array detector, is shown in Figure 3. In [26], it was shown analytically that the electric field of a wave at a distance z behind the axicon can be written in the paraxial approximation as Hereinafter, we use the coordinate system shown in Figure 2. We employed the polar coordinates (r , ϕ) for the axicon plane and the coordinates (ρ, θ) in the space behind the axicon and in the lens plane. The upper limit of the integral is the axicon radius a.
As shown above, the region of existence of quasi-Bessel beams is determined by expression (11). A simplified analytical expression for the electric field of the Bessel beam was derived from (13) for the region a << z < z A (z A = ka/κ): and for the intensity: One can see that the cross-section of the beam intensity exactly corresponds to the Bessel function. The analytical calculations in [19] were finished on this expression, and further results were obtained using numerical simulations and experiments. Since in this work, we are interested in the excitation of plasmons by means of diffraction of radiation on an end face of a cylindrical conductor, it is useful to derive, as an intermediate stage, analytical expressions for beams illuminating such a conductor. In the next section, we will find analytical expressions for converting beams described by expression (13) using a lens. Since the lens diameter is approximately equal to the axicon radius a, which is much larger than the diameter of the first rings of the Bessel beam, in the further calculations, we used expression (13), which does not contain further simplifications, as an initial expression for the wave illuminating the lens.

Transformation of Bessel Beams Formed by Binary Axicons into Perfect Vortex Beams Using a Lens
Below is the expression for a beam that has passed through a thin lens with a focal distance f [16,32]: where E(z L , ρ, θ) is the amplitude of the wavefield incident on the lens and a L is the lens radius. Placing a thin lens at the distance z = z L from the axicon and taking into account the fact that k − κ 2 /2k ≈ k z , after some transformations, we get the following expression at the distance z behind the lens: Here, r and ψ are the polar coordinates in the plane of the lens and a L is the radius of the lens. Since we have Making a trivial substitution ρ 2 = x in the integral with respect to ρ, introducing the variable Z L f = (1/z L + 1/z − 1/ f ) −1 , and expanding the upper limit of integration a L to infinity, we get the following integral, the value of which is given in [33] (see Appendix A): Substituting the resulting expression (20) into (19), we obtain the final expression for the field behind the lens: For the distance z = f, which is the main interest to us, we get the field amplitude where we omitted the insignificant constant phase factors in the exponent. The integrals in expressions (21) and (22) cannot be taken analytically; therefore, further calculations (see below) will be performed numerically. Based on the law of conservation of angular momentum, one can assume that the beam behind the lens has the same topological charge as the quasi-Bessel beam created by the axicon. In expressions (21) and (22), this dependence is hidden in the harmonic term under the integral (compare with Equation (5)).

Perfect Vortex Beams: Graphical Representation
As already mentioned in Section 3.1, Bessel beams formed by diffractive axicons have features that are not manifested in the intensity distribution, but are present in the phase distribution. When the beams are transformed by a lens, which is a Fourier transformer, these features are manifested in the spatial distribution too. Since our next task will be excitation of surface plasmon polaritons using perfect beams, the efficiency of plasmon excitation will depend on the features of the amplitude-phase distribution in the beams. In Figure 4, we present graphs of the distributions of fields and intensities in the lateral and longitudinal sections of perfect beams with the topological charges = 0, 1, 2 constructed based on Equations (21) and (22). Together with the plots for the Bessel beams presented in [26], they can serve as a useful reference for planning experiments with binary phase axicons. It is seen from the figure that in the case of beams created by binary spiral axicons, the intensity distribution at the focus is much more complicated than that for the ideal "perfect beams" (12). Instead of a single ring, we see (Figure 4a,b) a structure of nested rings (at = 0) or nested segments of spirals (at | | > 0). However, most energy is concentrated in the annular region corresponding to the radius R f . In Figure 4a, the direction of the twisting of the spirals corresponds to beams with a negative topological charge when the observer looks along the axis z. Figure 4c shows the effect of the diffraction of beams because of the limited beam diameter. It is seen that the beam cross-section is not symmetric with respect to the focal plane. This feature can be used for adjustment of the diameter of the "perfect beam" to the diameter of the conductor on which surface plasmon polaritons are to be excited. The phase distribution in the ring is also quite complex and can affect the efficiency of SPP generation.  I(r, 0, z) for 80 ≤ z ≤ 120 mm, calculated using Formula (21). The radiation wavelength is 141 µm.

Characteristics of SPPs on Axisymmetric Conductors
In this section, we will discuss the features of plasmon excitation by the end-fire coupling technique [13], which we must take into account when evaluating the applicability of one or another axicon for this. Surface plasmon polaritons [34] are essentially coupled oscillations of electrons in the surface layer of a conductor and a p-polarized electromagnetic field in an adjacent dielectric, which propagate along the surface. The field strength goes down on both sides of the interface, and the wave energy decreases due to ohmic losses in the conductor. In the case of cylindrical geometry, the normal component of the plasmon field value F(r) decreases according to the law (see, e.g., [35]) where I 1 (•) and K 1 (•) are modified first-order Bessel functions of the first and second kind, respectively, p m,d = k 2 s − k 2 0 ε m,d and ε m,d are the dielectric permittivities of the metal and dielectric, respectively, k s is the wave vector of plasmon, and a w is the radius of the cylinder.
Our prime interest is the efficiency of the excitation of the plasmon that rotates while propagating along a conductor with a frequency determined by the orbital angular momentum of the exciting beam; therefore, as a model of a plasmon, we will choose an electromagnetic wave that has the following field distribution over the surface and carries an orbital moment corresponding to the moment of the exciting wave. Here, γ is the decay constant of the amplitude of the surface wave during the propagation of the plasmon along the cylinder. In the terahertz range, for metals such as copper and gold, the plasmon propagation length varies from centimeters to tens of centimeters [34], that is, by an order of magnitude; the attenuation coefficient of the surface wave along the z-axis is 10 −2 mm −1 . With a conductor diameter of several millimeters, the electric field decays almost exponentially with respect to the radius: In this case, the decay constant β must be close to the respective value for plasmons on flat surfaces. For a gold surface with a zinc sulfide coating that is 0 to 1.5 µm thick, the measured β values are in the range of 0.3 ÷ 5 mm -1 [36]. According to experiments [37], the optimal thickness of the dielectric coating, which ensures good coupling of the plasmon with the surface and reduces parasitic radiation losses, is approximately 1 µm. In this case, the characteristic decay length for the evanescent wave in air is approximately 200-300 µm, that is, practically all energy of the plasmon field is concentrated within a cylindrical layer that is less than 0.6 mm thick.

Wavefields of Perfect Beams Obtained with Binary Axicons and Overlap Integral
In waveguide systems, the efficiency of the transformation of the exciting wavefield into a waveguide mode is usually estimated by the value of the overlap integral of the incident wavefield and the waveguide mode field (see, for example, [38]). The closer the distribution of the incident wave amplitude and phase to the expected distribution of these quantities in the waveguide mode, the greater the efficiency of excitation of the latter. Let us write the overlap integral for our case. It is known that the electric field of a plasmon on a cylinder has both radial and longitudinal components, while the incident wave at the lens focus has a plane wavefront; therefore, it is more convenient to write the overlap integral through the vectors of the magnetic field, which, in both cases, has only azimuthal components: Unlike the notation conventional for the bulk waves, as, for example, in [39], we write here not the total efficiency of the transformation of a free wave into a plasmon mode, but its derivative as a function of the azimuthal angle. This notation takes into account the difference between the excitation of a waveguide mode, for example, in a resonator, and the excitation of a plasmon on a surface. In the first case, the excited mode is spatial, and its distribution across the entire volume is known in advance. All the energy absorbed transmits into this mode. In the case of plasmons, the situation is different. A plasmon is excited at the conductor edge and moves in a certain direction. Simultaneously, another plasmon moving in parallel can be excited nearby with the same field distribution, but, generally speaking, with a different intensity, set by the local pumping intensity. In other words, in the first approximation, the excitation of plasmons in different regions can be considered independent. This is evidenced by the observation of parallel non-expanding tracks of plasmons created by the diffraction of a Bessel beam at the edge of a convex surface [40].
To calculate expression (27), one should substitute the expected distribution of the plasmon field in the air, and the field of the vortex beam illuminating the end of the cylinder should be substituted into it. The plasmon field is described by the expression F( f , r, ψ) = e −β(r−a w ) e i ψ , where the β value can be taken from the experimental data for the Au-ZnS-air surface [34], and the incident wavefield can be calculated according to (22). Figure 5 shows the intensity and phase distributions in the Bessel beams generated by binary axicons and in perfect beams in the focal plane of the lens. The numerical calculations were carried out using the diffraction integral [41], whereas the analytical calculations were done using Equation (22). The results of the calculations are very similar. The parameters of the axicons were selected in a way such that the diameter of the beams at the focus of the lens was approximately 7 mm, which is close to the diameter of the cylindrical conductors prepared for experiments. One can see that the width of the brightest part of the ring is about 0.6 mm, which is in good agreement with the plasmon decay length estimated above. The wavefields of perfect beams produced by binary axicons of different orders (see Figure 5) have variations of brightness over the azimuthal angle, which has a number equal to 2 . Thus, we can expect the generation of 2 parallel SPPs on a cylinder if the topological charge of a perfect beam is equal to .

Comparison of Perfect Beams Generated with Different Axicons
Since the field distributions in "perfect beams" created by binary axicons are quite far from ideal rings, we calculated the fields of Bessel beams and their images in the focal plane of a lens for all three axicons within the scalar diffraction theory. The calculation results are shown in Figure 6.
As seen from the third row in Figure 6, the calculated value of the ring radius in the focal plane of a lens with f = 100 mm for an ideal Bessel beam with a transverse wavenumber (4)) is 6.9 mm. The inner radius of all three rings in the figure is 6.6 mm. The radius of the maximum intensity of the ring created using the holo-

Comparison of Perfect Beams Generated with Different Axicons
Since the field distributions in "perfect beams" created by binary axicons are quite far from ideal rings, we calculated the fields of Bessel beams and their images in the focal plane of a lens for all three axicons within the scalar diffraction theory. The calculation results are shown in Figure 6.   As seen from the third row in Figure 6, the calculated value of the ring radius in the focal plane of a lens with f = 100 mm for an ideal Bessel beam with a transverse wavenumber κ = 3.1 (see Equation (4)) is 6.9 mm. The inner radius of all three rings in the figure is 6.6 mm. The radius of the maximum intensity of the ring created using the holographic axicon coincides well with the calculated value of 6.9 mm. The rest of the rings are much wider, and the numbers of spiral segments they have are 9 and 18, respectively. Using these results and adding to them the results of calculations for the same axicons with a topological charge of 3, we collected in Figure 7 the cross-sections of perfect beams produced with different axicons. Now, we summarize in Table 1 the main characteristics of perfect beams created using different types of axicons and discuss the results obtained.    = 100 mm). For other parameters, see Figure 6.

Discussion
The results of the analytical and numerical calculations of the characteristics of terahertz vortex beams created by the three types of axicons show that their characteristics depend substantially on the axicon type. Although all these axicons create quasi-Bessel beams, in which the amplitude-phase distributions of the field in the first ring are almost identical, the distributions of the field at the periphery differ-insignificantly, at first glance. However, when a lens focuses these beams, the main contribution into the formation of perfect beams, which are essentially the Fourier transform of a quasi-Bessel beam, is made by the beam periphery, in which the energy significantly exceeds the energy contained in the first ring. As a result, as seen in Figures 5-7 and Table 1, the shape of the annular PVB is highly dependent on the phase distribution profile of each axicon. The diffraction phenomena associated with the relatively small number of zones in the axicons, which should be rather wide in the terahertz range because of the large wavelength, may also play a certain role in the formation of the PVB profile.
From this, we can draw two conclusions. First, perfect beams formed in different experiments always bear the imprint of their origin, and the properties of each of them are largely individual. The second concerns the Bessel beams themselves. It is very likely that their inherent self-recovery property after passing through inhomogeneities also depends on their origin, since self-recovery is provided by the arrival of converging plane waves from the beam periphery into the optical axis region. Apparently, this circumstance should be taken into account in propagation of Bessel beams in inhomogeneous media.
Another issue to discuss is the dependence of PVB parameters on the magnitude of the topological charge. In [27], it was shown theoretically and experimentally that, in contrast to the conclusions of the simplified theory, the radius and width of the rings of PVBs created in the visible range using a spatial light modulator grow with the topological charge. The results of our calculations show that in the case of terahertz vortices created using diffraction axicons, we observe no strong dependence on the topological charge, although there is a slight upward trend. We can trust the accuracy of calculations by this program, since, in [26], it was demonstrated that the images of perfect beams obtained in the experiment coincided with the calculated ones.
In conclusion, let us discuss the possibilities of using axicons to excite surface plasmon polaritons on cylindrical waveguides. A factor favorable for the excitation of terahertz plasmons is the very small change in the radius and width of the rings of perfect plasmons with a change in the topological charge. Apparently, radially polarized PVBs can be applied to the excitation of azimuth-uniform SPPs using a holographic axicon. In this case, the half-width of the PVB ring is approximately 400 µm, which is in good agreement with the decay length of plasmons on a gold surface covered with a zinc sulfide layer that is about 0.5 µm thick. When using other axicons, one can expect the separation of SPPs in azimuth into or 2 plasmons moving in parallel, which, however, carry the same total orbital angular momentum. Since in this case, the half-width of the rings is of the order of 1 mm, they can be used at a dielectric coating thickness of about 0.1 µm, at which the decay length is close to this value. With linear polarization of the PVB, plasmons will be excited on the opposite sides of the cylinder, where the wave electric field vector is normal to the surface. Finally, note that the limited aperture of the beam on the axicon, which causes diffraction of the beam, can also be beneficial in experiments. Figure 4 shows that near the focal plane, the ring diameter changes smoothly; by moving the end face of the cylinder along the axis, one can optimize the efficiency of the plasmon excitation due to change in the overlap integral value (27).