## Abstract

We present the generation of arbitrary order Bessel beams at 0.3 THz through the implementation of suitably designed axicons based on 3D printing technology. The helical axicons, which possess thickness gradients in both radial and azimuthal directions, can convert the incident Gaussian beam into a high-order Bessel beam with spiral phase structure. The evolution of the generated Bessel beams are characterized experimentally with a three-dimensional field scanner. Moreover, the topological charges carried by the high-order Bessel beams are determined by the fork-like interferograms. This 3D-printing-based Bessel beam generation technique is useful not only for THz imaging systems with zero-order Bessel beams but also for future orbital-angular-momentum-based THz free-space communication with higher-order Bessel beams.

© 2015 Optical Society of America

## 1. INTRODUCTION

Terahertz (THz) waves have exhibited strong potential for applications in the field of imaging [1] for the sake of its penetration capacity for many optically opaque materials, such as wood, plastic, and biological tissue. Conventional THz imaging techniques usually employ dielectric convex lenses or off-axis parabolic mirrors to focus the incident THz beam. According to Rayleigh’s criteria, the fixed length of the focusing elements will produce an inherent trade-off between the lateral resolution (i.e., focal spot) and axial resolution (i.e., focal depth). However, in many potential applications, such as stand-off personal screening or defects detection with the uncertain location of thick objects, a tightly focused THz beam that retains its size over a large depth of focus is highly desirable. Recently, Bessel beams have attracted much attention for their potential to overcome this trade-off limitation [2,3]. Through the use of a simple conical-shaped lens known as an axicon, Bessel beams can be generated efficiently [4]. By adjusting the apex angle of the axicon and the diameter of the incident Gaussian beam, one can get a Bessel beam ensuring continuous highly invariant lateral resolution over an extended axial length, which is attractive for THz imaging applications. In recent years, intensive research efforts have been performed conforming the improvement of Bessel beams on the depth of field in two-dimensional THz imaging systems using continuous THz waves [5 –7], as well as pulsed THz radiation [8,9]. Moreover, due to their propagation-invariant lines of focus over an extended depth-of-field, Bessel beams have the potential to fulfill the condition of straight ray lines of the Radon transform in THz-computed tomography, yielding a considerable improvement in the 3D reconstruction image of an object [10].

In practice, ideal Bessel beams cannot be created due to their infinite lateral extent and infinite amount of energy, while quasi-Bessel beams which possess the nondiffraction properties over a finite range can be experimentally realized with the classical axicons, proposed by McLeod in 1954 [11]. In our work, the experimentally generated quasi-Bessel beams will be referred to simply as Bessel beams. Besides the regular bulk axicon which is based on pure refraction, novel types of axicons using the Fresnel-type approach and holographic methods have also been proposed. By applying the Fresnel lens concept to a refractive axicon, an equivalent of the refractive axicon can be achieved, called a Fresnel axicon (fraxicon) [12]. Holographic optical elements can mimic the phase distribution of the refractive axicon, with on-axis configuration or off-axis configuration [13 –17]. It is worthy to note, with programmable spatial light modulators (SLMs), such as digital micromirror (DMD) [18] and liquid crystal (LC) systems [15], Bessel beam generation can be dynamically controlled in the visible and infrared wavebands. However, in the terahertz regime, those commercial SLMs do not operate because of the lack of materials with the desired terahertz modulation. In the microwave and THz bands, some approaches like binary hologram [19,20], holographic metasurface [21], and metamaterial lens [22] have been demonstrated to produce Bessel beams. So far, the generation of THz Bessel beams has mostly relied on conventional bulk refractive axicons, which are manufactured through the use of traditional computer numerical control (CNC) milling techniques [7,23,24]. Currently, there is a growing interest in using 3D printing technology to make THz devices, such as waveguides [25], electromagnetic bandgap structures [26], and computer-generated volume holograms [27]. In our earlier work, we reported the generation of THz vortex beams carrying orbital angular momentum (OAM) states using 3D printed spiral phase plates (SPPs) [28,29], which is shown to be more straightforward and efficient compared with the mechanically polishing method reported previously [30 –32]. So far, to the best of our knowledge, studies via 3D printed elements to generate THz Bessel beams have not been reported.

In this paper, we extend the 3D printing technique to fabricate refractive elements for generating arbitrary order Bessel beams at the THz frequency. Special attention has been paid to minimize the volume of the elements. By changing the base angle of the fabricated axicons, the central spot size of the zero-order Bessel beams are modulated. Thus, the axicon-generated zero-order Bessel beams can meet the specific demands in practical THz imaging situations. Moreover, helical axicons are designed to generate higher-order Bessel beams, which can carry OAM states over extended distances in a nondiffracting manner. Our proposed Bessel beam generation technique may be useful in the field of THz imaging applications as well as OAM-based THz free-space optical (FSO) communication [33,34].

## 2. DESIGN AND FABRICATION OF THE PHASE PLATES

In 1987, Durnin first proposed the Bessel beam as a special nondiffracting solution of the wave equation in free space [2]. By solving the Helmholtz equation in the circular cylindrical coordinate system, the electric field of an $l$-order Bessel beam is given by

where ${J}_{l}$ is the $l$-order Bessel function of the first kind; $\rho $ and $\phi $ are the radial and azimuthal coordinates in the transverse plane, respectively; ${k}_{\rho}$ and ${k}_{z}$ are the radial and longitudinal components of the free-space wavevector, respectively, such that ${k}_{\rho}^{2}+{k}_{z}^{2}={k}^{2}$. According to Eq. (1), the lateral electric field distribution is described in terms of the $l$-order Bessel function of the first kind. Due to the spatial phase distribution $\mathrm{exp}(-jl\phi )$, the $l$-order Bessel beams carry an OAM state with the topological charge of $l$.In general, OAM beams can be produced with a spiral phase plate (SPP) illuminated by the Gaussian beam, as shown schematically in Fig. 1(a). A typical SPP is composed of a medium of refractive index $n$ with an azimuthally varying thickness, given by ${h}_{1}(\phi )=l\lambda \phi /(n-1)2\pi $. With one surface set flat, the required thickness of the SPP can be obtained simply by controlling the height of the other surface. The total height of the SPP is chosen as $\mathrm{\Delta}h=l\lambda /(n-1)$, so that the total phase delay around the plate center is an integer multiple of $2\pi $, i.e., $2\pi l$. In analogy with the SPP which induces an azimuthal $l\phi $ phase ramp, the axicon can be called as the conical phase plate, since it introduces a radial phase delay. As shown in Fig. 1(b), the axicon can convert an incident Gaussian beam into a zero-order Bessel beam. The thickness of a typical axicon keeps constant in the azimuthal direction and increases in proportion to the radius toward the center, given by ${h}_{2}(\rho )=\rho \text{\hspace{0.17em}}\mathrm{tan}(\gamma )$. By replacing the input Gaussian beam with an OAM beam, which has a helical wavefront, one can also use the axicon to produce higher-order Bessel beams [35]. Moreover, by simply combining the axicon with the SPP, a new element called a helical axicon can be achieved, which has double height-varying surfaces. As shown in Fig. 1(c), the helical axicon can directly transform the Gaussian beams into high-order Bessel beams. A similar concept has been mentioned in the optical spectral range [36,37].

Figure 2(a) shows the schematic of the SPP designed to produce an OAM beam with $l=1$. Both the magnitude and sign of the value $l$ can be controlled by increasing the total height $\mathrm{\Delta}h$ of the SPP in the clockwise or counterclockwise direction. It is obvious that the larger the value $l$ is, the thicker the SPP is, and the heavier weight, as well as greater loss will arise. Utilizing the $2\pi $ phase jump, we remove the surplus thickness corresponding to multiple $2\pi $ phase changes and thus reduce the thickness of the $l>1$ SPP to the wavelength scale. The SPP with topological charge $l$ contains $l$ equal modules which have an identical thickness $\mathrm{\Delta}h=\lambda /(n-1)$. As a demonstration, Fig. 2(b) shows the scheme of an $l=6$ SPP. Six basic $l=1$ SPPs are compressed from an angular range of [0, $2\pi $] to [0, $\pi /3$] and pushed together to constitute the $l=6$ SPP. This idea has been applied in our previous work for designing high-order SPPs and can be further extended to the helical axicons in this work. Figure 2(c) gives the schematic drawing of a typical axicon, whose height is determined by the base diameter and base angle, as depicted in Fig. 1(b). By combing the SPP and axicon, the helical axicon can be created, as shown in Fig. 1(c). Furthermore, much more compact helical axicons can be achieved by setting one surface flat and changing the height of the other surface in azimuthal and radial directions simultaneously. The scheme of the $l=6$ helical axicon is shown in Fig. 2(d). To differentiate the geometry of the helical axicon and the SPP with the same topological charge $l$, we pay attention to the cross section of each module of the two phase plates. Shown as the red regions in Figs. 2(b) and 2(d), the SPP has a cross section of a rectangle, while the helical axicon has a cross section of a right triangle. Each module of both the SPP and the helical axicon increases its height azimuthally to the angle of $2\pi /l$, thus imprinting a total phase shift of $2\pi $ on the incident wave. The thickness of each module of the helical axicon is given by

For the special case $l=0$, the helical axicon evolves to be an axicon. Applying small angle approximations, $\mathrm{tan}(\gamma )=\gamma $ and ${k}_{\rho}=k(n-1)\gamma $, the transmittance functions of the SPP, axicon, and helical axicon, respectively, are given by

Obviously, the transmittance function of the helical axicon is the product of that of the SPP and axicon, namely, ${T}_{3}(\phi ,\rho )={T}_{1}(\phi ,\rho {)}^{*}{T}_{2}(\phi ,\rho )$. It can be derived that the phase modulation induced by the helical axicon is the sum of the SPP and axicon, in the form

With ray optics analysis, the first term of the above formula represents the radial part which depends linearly on the radius and results in conical refraction. The second term represents the azimuthal part which depends linearly on the azimuthal angle and results in azimuthal refraction. Accordingly, the base angle $\gamma $ together with the order $l$ specify a helical axicon. In our work, we fabricated six axicons with $\gamma $ ranging from 5° to 30° (shown in the first row of Fig. 3) and six helical axicons with $l$ ranging from 1 to 6 (shown in the second row of Fig. 3). All the helical axicons have a fixed base angle $\gamma =10\xb0$ to make a consistent comparison with the $\gamma =10\xb0$ axicon. Six SPPs with $l$ ranging from 1 to 6 were also fabricated (shown in the third row of Fig. 3) for the comparison with the helical axicons in the following discussions.

All the phase plates were manufactured by employing a commercial Objet30 3D printer with the composition polymer ink “FullCure835VeroWhite” material. Before the fabrication process, the optical properties of the polymer material have been characterized with a Zomega-Z3 THz time-domain spectrometer (THz-TDS). At the frequency of 300 GHz, the refractive index and the absorption coefficient of the polymer material are about 1.655 and $1.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, respectively. All the phase plates have a diameter of 52.8 mm. As can be deduced from ${h}_{2}(\rho )=\rho \text{\hspace{0.17em}}\mathrm{tan}(\gamma )$, the maximum thickness changes of the six axicons are 2.31, 4.66, 7.07, 9.61, 12.31, and 15.24 mm, respectively. According to $\mathrm{\Delta}h=\lambda /(n-1)$, all the SPPs have an equal maximum thickness of 1.53 mm. The maximum thickness of helical axicons is 6.19 mm, which is equivalent to that of a $\gamma =10\xb0$ axicon and the SPP combined. For mounting the phase plates on the holder, we designed an additional base for every phase plate with a diameter of 56.8 mm and a thickness of 2.2 mm. It is worth remarking that the surface smoothness of the phase plates is determined by the printing precision. The Objet30 3D printer has a resolution of 600 dpi (42 μm) in the $x$ and $y$ directions and 900 dpi (28 μm) in the $z$ direction.

## 3. BESSEL BEAM PARAMETERS AND MEASUREMENT SYSTEM

As illustrated in Fig. 1(b), the light from the incident Gaussian beam refracts toward the central $z$ axis after passing through the axicon, and its interference defines a Bessel beam. The properties of the axicon-produced Bessel beam are determined by the base angle $\gamma $ of the axicon and the beam waist radius ${\omega}_{0}$ of the incident Gaussian beam. Using Snell’s law, the semi-apex angle ${\alpha}_{0}$ by which the incoming ray is refracted, can be calculated as

Here, the apex angle of the axicon is defined as $\tau =\pi -2\gamma $, and $n$ and ${n}_{0}$ are the refractive index of the axicon and the surrounding medium, respectively. In our work, the ambient medium is air (${n}_{0}=1$). For a small base angle $\gamma $, the above expression can be simplified to obtain the approximation formula as

Such a beam propagates diffraction-free over a distance ${z}_{\mathrm{max}}$, starting from the tip of the axicon. As shown in the geometrical configuration of Fig. 1, ${z}_{\mathrm{max}}$ can be obtained from Snell’s law and is estimated as

Figure 4 shows schematically the experimental setup for obtaining and measuring the generated Bessel beams at THz frequencies. A 100 GHz Gunn diode (Spacek Labs model GW-102P) coupled to a frequency tripler delivers 0.3 mW at 300 GHz. Then the emitted THz beam is collimated by a TPX lens and directed onto the phase plates. In the schematic diagram below, an axicon is given as an example of the phase plates. The inset in the upper right corner of Fig. 4 illustrates the 2D spatial profile of the collimated THz beam in the input plane of the phase plates. The intensity of the incident beam shows a Gaussian-shaped profile with a diameter of $2{\omega}_{0}=46\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ measured at full width at half-maximum (FWHM). The diameter of the collimated THz beam is smaller than that of the phase plates, which ensures the high conversion efficiency. Absorbing material is placed around the holder to terminate the incident radiation outside of the phase plates. A Schottky diode (Virginia Diodes Inc.) mounted on three-dimensional motorized translational stages, is employed to measure the resulting spatial intensity profiles after the phase plates. Point-to-point scanning is performed by the motorized stages along $x$ and $y$ axes with a step resolution of 0.5 mm, over a transversal coverage area of $90\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}*90\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. The detector can be further moved along the $z$ axis with a maximum range of 300 mm. Before the measurement, we have checked the transmitter–detector alignment to ensure that the beam center is accordant with the center of the translation-stage-scanning range. Lock-in detection with a 300 Hz mechanical chopper placed after the transmitter is used to reduce noise and extract a reliable signal. The elements in the dashed box are used for obtaining the interference pattern, which will be discussed in detail later. By moving away those elements, the spatial intensity distribution of the generated beams after the phase plates can be directly measured.

## 4. EXPERIMENT RESULTS AND DISCUSSION

#### A. Characterization of Zero-Order Bessel Beams

Utilizing the configuration in Fig. 4, we investigate the detailed features of the generated Bessel beams by switching the six axicons sequentially. The base of the axicons (the input plane) is located at $z=-16\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, and the center is at $x=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, $y=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. First, we measured the intensity distributions along the propagation direction in the $x\u2013z$ plane (longitudinal profiles) over a coverage area of $90\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}*300\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, as summarized in the first row of Fig. 5(c). One can clearly see that a beam is formed along the $z$ axis, which shows a central peak and periodic side-lobes. To investigate more features of the beams, the axial intensity distributions corresponding to different $\gamma $ are extracted and plotted together for comparison in Fig. 5(b). We define the diffraction-free range ${z}_{\mathrm{max}}$ as the interval from the apex tip to half the maximum intensity along the propagation direction. Then the corresponding values of ${z}_{\mathrm{max}}$ are plotted together as discrete triangles in Fig. 5(a). The theoretical values of ${z}_{\mathrm{max}}$ are calculated according to Eq. (7) using the parameter $2{\omega}_{0}=46\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ and plotted as the black curve in Fig. 5(a). Since the theoretical ${z}_{\mathrm{max}}$ value of the $\gamma =5\xb0$ axicon is beyond the maximum range of the motorized stage in the $z$ axis which has a value of 300 mm, the experimental ${z}_{\mathrm{max}}$ value of the $\gamma =5\xb0$ axicons is not presented in the figure. As can be seen, the experimentally obtained values of ${z}_{\mathrm{max}}$ match well to the theoretically calculated curve.

At the axial distances ${z}_{f}$ [marked as the white dotted lines in the first row of Fig. 5(c)], we measure the corresponding cross-sectional intensity distributions in the $x\u2013y$ plane (transverse profiles). The normalized transverse profiles are summarized in the second row of Fig. 5(c). The profiles show an axially symmetric intensity pattern, which consists of a bright core surrounded by concentric rings of decreasing radial intensity. Along with the increase of $\gamma $, the size of the central spot decreases and more rings arise. At the distance of ${z}_{f}$, the corresponding 1D cross-sectional profiles through the center of the $x\u2013y$ profile along the $x$ axis are plotted as indicated by the dotted red line in the third row of Fig. 5(c). The radial profiles are well fitted with the absolute square of the Bessel function of the first kind ${J}_{0}^{2}({k}_{\rho}^{*}\rho )$ (marked as the solid black line), and the corresponding Bessel functions are indicated in the legend. Since zero-order Bessel function ${J}_{0}({k}_{\rho}^{*}\rho )$ has its first zero at ${k}_{\rho}^{*}\rho =2.4048$, the radius of the central spot is given by $2.4048/{k}_{\rho}$. When ${k}_{\rho}=0$, the Bessel beam reduces to a plane wave, and when ${k}_{\rho}=k$, it reaches the minimum diameter of about ¾ $\lambda $. With regard to the measured results, the beam core size can be defined as the distance from the center of the beam to the first radial intensity minimum and is found to decrease with increasing $\gamma $.

It can be concluded that a larger base angle leads to a smaller spot of Bessel beam, but a much shorter nondiffraction zone. Therefore, by properly adjusting the base angle of the axicon, we can fine-tune the focal length and focal depth of the generated Bessel beam, to fulfill the specific requirements in practical terahertz imaging situations.

#### B. Comparison between High-Order Bessel Beams and OAM Beams

We further investigate the higher-order Bessel beams generated by helical axicons with $l$ ranging from 1 to 6. Figure 6(a) shows the $x\u2013z$ longitudinal profiles and the corresponding $x\u2013y$ transverse profiles. As will be discussed in detail later, all the longitudinal profiles of the high-order Bessel beams have maximum values around ${z}_{f}=132\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. Therefore, all the $x\u2013y$ transverse profiles are measured at the distance ${z}_{f}=132\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, consistent with that of the $\gamma =10\xb0$ axicon. We can clearly see from the transverse and longitudinal profiles that the high-order Bessel beam has a central dark region surrounded by bright rings, and this form remains invariant over a propagation distance.

For comparative purposes, we also measure the longitudinal and transverse profiles of the OAM beams generated by SPPs with $l$ ranging from 1 to 6, shown in the first and second row of Fig. 6(b), respectively. As can be seen, both the high-order Bessel beams and the OAM beams show a central dark line in $x\u2013z$ longitudinal profiles, corresponding to the dark core in the $x\u2013y$ transverse profiles. Furthermore, the diameter of the central dark zone, namely, the inter ring of the high-order Bessel beams or the dark core of the OAM beams, increases with the order $l$. In analogy with the zero-order Bessel beam, the high-order Bessel beams show properties of nondiffraction, where the size of the inner ring remained invariant over the propagation distance. In contrast, the OAM beams exhibit strong divergence with the dark core expanded along the $z$ axis. From the $x\u2013z$ longitudinal profiles of OAM beams with different values of $l$ [the first row of Fig. 6(b)], it can be found that the OAM beams with larger $l$ diverge more. According to Refs. [38,39], the beam divergence should scale linearly with the square root of $l$, since the beam waist of the Gaussian term is kept constant for different SPPs, in our case where the OAM beams are produced by the SPPs.

Then we pay our attention to the transverse profiles at a specific propagation distance. At the distance of ${z}_{f}=132\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, the corresponding 1D cross-sectional profiles along the $x$ axis are plotted in Figs. 7(a) and 7(b). The radial profiles agree well with the absolute square of the $l$-order Bessel function of the first kind ${J}_{l}^{2}({k}_{\rho}^{*}\rho )$ (the solid black line), and the corresponding Bessel functions are indicated in the legend. To further investigate the nondiffraction characteristics, we extract their maximal intensity values at each distance from the longitudinal profiles, and the results are presented in Figs. 7(c) and 7(d), respectively. The maximal intensity of high-order Bessel beams evolves in a manner analogous to that of the zero-order Bessel beam generated by the $\gamma =10\xb0$ axicon [the red line in Fig. 5(b)], which reaches the maximum values around the distance of ${z}_{f}=132\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. In contrast, the OAM beams show descended maximum values along the propagation direction, since the hollow size expands monotonically due to the diffraction effect.

The discussions above indicate that, under the same Gaussian beam illumination, the $l$-order Bessel beams generated by helical axicons show a much smaller size compared to the $l$-order OAM beams generated by SPPs, at the same propagation distance $z$. Furthermore, the higher-order Bessel beams remain invariant to the inner ring size over a propagation distance, indicating that they not only carry the OAM states but still have the propagation-invariance property. Hence, high-order Bessel beams can relax the stringent requirements on the size of aperture at the receiving optical system when used for communication purposes. It can be concluded that high-order Bessel beams are advantageous as OAM carriers in long-distance OAM-based FSO communication systems.

#### C. 3D Evolution of Bessel Beams and OAM Beams

To obtain deeper insight into the evolution properties, the volumetric renderings of the generated beams are investigated. We record each $x\u2013y$ transverse intensity profile while varying the $z$ axis distance over a range of 0 to 300 mm in steps of 20 mm. Hence, 3D intensity data are taken. With the use of 3D visualization software Voxler (Golden Software Inc.), we can conveniently reconstruct and visualize the beams from the 3D intensity data. Figure 8 shows the 3D visualization of the zero-order Bessel beam, first-order Bessel beam, and $l=1$ OAM beam. The first row of Fig. 8 gives the volume rending with the alpha blending rendering method, while the second row gives the transverse and longitudinal cross-sectional intensity profiles.

As shown in Figs. 8(a) and 8(b), it is obvious that the first-order Bessel beam behaves similarly to the zero-order Bessel beam along the propagation distance. Both of them have similar locations of maximum values and similar diffraction ranges. Furthermore, from the cross-sectional intensity profiles of the first-order Bessel beam in Fig. 8(e), a hollow-center intensity distribution can be distinguished. As a contrast, the zero-order Bessel beam has a central nonzero intensity, as can be seen from Fig. 8(d). Then we compare the first-order Bessel beam with the $l=1$ OAM beam. The OAM beam shows as a size-expanding hollow cylinder. In contrast, the high-order Bessel beam maintains the “diffraction-free” profile and keeps an invariant hollow center. Thus, they are advantageous as OAM carriers for signal delivery to distance receivers.

In the optical regime, the analytical expression of the intensity distribution behind an axicon can be obtained based on scalar diffraction theory [40], and the treatments ignore diffraction effects on the axicon edges. However, in our case of the THz regime, the phase plates have a diameter which is only one order of magnitude greater than the wavelength and slightly larger than the waist of the incident beam. The edge diffraction effects caused by truncation at the edge of the phase plates are pronounced and have an important influence on the intensity distribution. So the scalar diffraction analysis is not suitable here. A full vector field treatment of the electromagnetic field is needed for more precise description of the propagation properties of the generated beam. Numerical simulation based on a rigorous FDTD calculation can give a better explanation of the experimental results [5].

#### D. Measurement of the OAM States Carried by the High-Order Bessel Beams

For the application of high-order Bessel beams in future FSO communications, it is of great importance to efficiently detect the OAM information from the beams. The most common technique to observe the OAM state is based on the fork-shaped interference pattern, by overlapping the OAM beam with a tilted reference plane wave. The number of prongs of the interference fork indicates the magnitude of topological charge $l$. The sign of $l$ is implied from the prongs up or down direction. This interferometric measurement is realized via the experimental setup in Fig. 4. The silicon wafer in the dashed box of Fig. 4 serves as the beam splitter, which offers a reference beam. With the combination of the elements in the dashed box, the interference patterns of the reference beam with the Gaussian beam, zero-order Bessel beams are generated by the $\gamma =10\xb0$ axicon and high-order Bessel beams with $l=1$, 2, 3 are recorded and shown in Figs. 9(a)–9(c), respectively.

Compared with the Gaussian beam [Fig. 9(a)], the interference pattern of the zero-order Bessel beam [the first row of Fig. 9(b)] shows triangular stripes at the center. The interference stripes bend to the right, which is associated with the conical wavefront of the zero-order Bessel beam. As shown in the interference patterns of high-order Bessel beams [the first row of Fig. 9(c)], the interference stripes show a fork-like splitting, and the number of the splitting scales is $l+1$. The direction of the splitting (“fork up”) indicates that all the high-order Bessel beams have the same direction of the azimuthal phase slope. One can also notice that the interference stripes bend the same amount as that in Fig. 9(b), which is determined by the base angle $\gamma $ of the helical axicon. If we reverse the helical axicons with the convex surface toward the incident beam, the interference patterns [shown in the second row of Fig. 9(c)] display reverse fork directions (“fork down”), implying that the topological charges of the generated Bessel beams change sign. However, the bending direction and extent of the interference stripes remain unchanged, indicating the radial phase slopes of the beams do not change. Those interferograms are also numerically calculated using the phase distribution of Eq. (4), shown in the third row of Fig. 9, which agree well with the measured results (the second row of Fig. 9). Our experimental results confirm the OAM states carried by the high-order Bessel beams generated by the helical axicons.

More recently, novel proposals have been reported for the detection of OAM-carrying high-order Bessel beams by using the digital axicon [41] and the mode sorter [42], both of which are programmed on SLMs in the optical regime. It is worth pointing out that, in millimeter wave and THz bands, those phase patterns can be conveniently realized with freeform refractive optical elements by 3D fabrication [39,43].

## 5. CONCLUSION

In this article, we have demonstrated the Bessel beam generation based on 3D printed phase plates at 0.3 THz. The generated Bessel beams are experimentally characterized by the transverse and longitudinal cross-sectional intensity profiles, as well as the 3D volume-rendered intensity distributions. Moreover, the OAM states carried by high-order Bessel beams are measured by the fork-like interferograms. Our proposed 3D printing method makes the realization of Bessel beams reliable and flexible, which is appropriate for the THz imaging system with zero-order Bessel beams and the THz FSO communication system with higher-order Bessel beams.

Since the polymer material exhibits high transmission, the 3D printed bulk devices can serve as suitable quasi-optical components for wavefront shaping in the THz regime. This concept can be extended to generate other types of THz exotic beams, such as Airy beams associated with continuous cubic phase distribution [44,45].

## Funding

National Natural Science Foundation of China (NSFC) (11574105, 61177095, 61405063, 61475054); Natural Science Foundation of Hubei Province (2012FFA074, 2013BAA002); Shanghai Academy of Spaceflight Technology (SAST) (201345); Technology Innovation Foundation From Innovation Institute of Huazhong University of Science and Technology (CX14-070, CXY13M009, CXY13Q015); the Fundamental Research Funds for the Central Universities (2014QN023, 2014ZZGH021); Wuhan applied basic research project (20140101010009).

## REFERENCES

**1. **M. Tonouchi, “Cutting-edge terahertz technology,” Nat. Photonics **1**, 97–105 (2007). [CrossRef]

**2. **J. Durnin, “Exact solutions for nondiffracting beams. The scalar theory,” J. Opt. Soc. Am. A **4**, 651–654 (1987). [CrossRef]

**3. **D. McGloin and K. Dholakia, “Bessel beams: diffraction in a new light,” Contemp. Phys. **46**, 15–28 (2005). [CrossRef]

**4. **S. Monk, J. Arlt, D. A. Robertson, J. Courtial, and M. J. Padgett, “The generation of Bessel beams at millimetre-wave frequencies by use of an axicon,” Opt. Commun. **170**, 213–215 (1999). [CrossRef]

**5. **G. Ok, S.-W. Choi, K. H. Park, and H. S. Chun, “Foreign object detection by sub-terahertz quasi-Bessel beam imaging,” Sensors **13**, 71–85 (2012). [CrossRef]

**6. **C. Liu, X. Wei, Z. Zhang, K. Wang, Z. Yang, and J. Liu, “Terahertz imaging system based on Bessel beams via 3D printed axicons at 100 GHz,” Proc. SPIE **9275**, 92751Q (2014).

**7. **Z. Zhang and T. Buma, “Terahertz imaging in dielectric media with quasi-Bessel beams,” Proc. SPIE **7938**, 793806 (2011). [CrossRef]

**8. **A. Bitman, I. Moshe, and Z. Zalevsky, “Improving depth-of field in broadband THz beams using nondiffractive Bessel beams,” Opt. Lett. **37**, 4164–4166 (2012). [CrossRef]

**9. **A. Bitman, I. Moshe, and Z. Zalevsky, “Broadband THz, extended depth of focus imaging based on step phase mask aided interferometry,” Opt. Commun. **309**, 1–5 (2013). [CrossRef]

**10. **A. Bitman, S. Goldring, I. Moshe, and Z. Zalevsky, “Computed tomography using broadband Bessel THz beams and phase contrast,” Opt. Lett. **39**, 1925–1928 (2014). [CrossRef]

**11. **J. H. McLeod, “The axicon: a new type of optical element,” J. Opt. Soc. Am. **44**, 592 (1954). [CrossRef]

**12. **I. Golub, “Fresnel axicon,” Opt. Lett. **31**, 1890–1892 (2006). [CrossRef]

**13. **J. Turunen, A. Vasara, and A. T. Friberg, “Holographic generation of diffraction-free beams,” Appl. Opt. **27**, 3959–3962 (1988). [CrossRef]

**14. **A. Vasara, J. Turunen, and A. T. Friberg, “Realization of general nondiffracting beams with computer-generated holograms,” J. Opt. Soc. Am. A **6**, 1748–1754 (1989). [CrossRef]

**15. **J. A. Davis, J. Guertin, and D. M. Cottrell, “Diffraction-free beams generated with programmable spatial light modulators,” Appl. Opt. **32**, 6368–6370 (1993). [CrossRef]

**16. **J. Leach, G. M. Gibson, M. J. Padgett, E. Esposito, G. McConnell, A. J. Wright, and J. M. Girkin, “Generation of achromatic Bessel beams using a compensated spatial light modulator,” Opt. Express **14**, 5581–5587 (2006). [CrossRef]

**17. **M. Bock, S. Das, and R. Grunwald, “Programmable ultrashort-pulsed flying images,” Opt. Express **17**, 7465–7478 (2009). [CrossRef]

**18. **L. Gong, Y.-X. Ren, G.-S. Xue, Q.-C. Wang, J.-H. Zhou, M.-C. Zhong, Z.-Q. Wang, and Y.-M. Li, “Generation of nondiffracting Bessel beam using digital micromirror device,” Appl. Opt. **52**, 4566–4575 (2013). [CrossRef]

**19. **J. Salo, J. Meltaus, E. Noponen, J. Westerholm, M. Salomaa, A. Lönnqvist, J. Säily, J. Häkli, J. Ala-Laurinaho, and A. Räisänen, “Millimetre-wave Bessel beams using computer holograms,” Electron. Lett. **37**, 834–835 (2001). [CrossRef]

**20. **J. Meltaus, J. Salo, E. Noponen, M. M. Salomaa, V. Viikari, A. Lönnqvist, T. Koskinen, J. Säily, J. Häkli, and J. Ala-Laurinaho, “Millimeter-wave beam shaping using holograms,” IEEE Trans. Microwave Theor. Tech. **51**, 1274–1280 (2003). [CrossRef]

**21. **B. G. Cai, Y. B. Li, W. X. Jiang, Q. Cheng, and T. J. Cui, “Generation of spatial Bessel beams using holographic metasurface,” Opt. Express **23**, 7593–7601 (2015). [CrossRef]

**22. **M. Q. Qi, W. X. Tang, and T. J. Cui, “A broadband Bessel beam launcher using metamaterial lens,” Sci. Rep. **5**, 11732 (2015). [CrossRef]

**23. **N. Trappe, R. Mahon, W. Lanigan, J. A. Murphy, and S. Withington, “The quasi-optical analysis of Bessel beams in the far infrared,” Infrared Phys. Technol. **46**, 233–247 (2005). [CrossRef]

**24. **M. Shaukat, P. Dean, S. Khanna, M. Lachab, S. Chakraborty, E. Linfield, and A. Davies, “Generation of Bessel beams using a terahertz quantum cascade laser,” Opt. Lett. **34**, 1030–1032 (2009). [CrossRef]

**25. **S. Pandey, B. Gupta, and A. Nahata, “Terahertz plasmonic waveguides created via 3D printing,” Opt. Express **21**, 24422–24430 (2013). [CrossRef]

**26. **Z. Wu, J. Kinast, M. Gehm, and H. Xin, “Rapid and inexpensive fabrication of terahertz electromagnetic bandgap structures,” Opt. Express **16**, 16442–16451 (2008). [CrossRef]

**27. **W.-R. Ng, D. R. Golish, H. Xin, and M. E. Gehm, “Direct rapid-prototyping fabrication of computer-generated volume holograms in the millimeter-wave and terahertz regime,” Opt. Express **22**, 3349–3355 (2014). [CrossRef]

**28. **X. Wei, “Orbit angular momentum multiplexing in 0.1-THz free-space communication via 3D printed spiral phase plates,” in *Proceedings of Conference on Lasers and Electro-Optics (CLEO)* (2014), paper STu2F.2.

**29. **L. Zhu, X. Wei, J. Wang, Z. Zhang, Z. Li, H. Zhang, S. Li, K. Wang, and J. Liu, “Experimental demonstration of basic functionalities for 0.1-THz orbital angular momentum (OAM) communications,” in *Proceedings of Optical Fiber Communication Conference (OFC)* (2014), paper M3K.7.

**30. **P. Schemmel, G. Pisano, and B. Maffei, “A modular spiral phase plate design for orbital angular momentum generation at millimetre wavelengths,” Opt. Express **22**, 14712–14726 (2014). [CrossRef]

**31. **P. Schemmel, S. Maccalli, G. Pisano, B. Maffei, and M. W. R. Ng, “Three-dimensional measurements of a millimeter wave orbital angular momentum vortex,” Opt. Lett. **39**, 626–629 (2014). [CrossRef]

**32. **K. Miyamoto, K. Suizu, T. Akiba, and T. Omatsu, “Direct observation of the topological charge of a terahertz vortex beam generated by a Tsurupica spiral phase plate,” Appl. Phys. Lett. **104**, 261104 (2014). [CrossRef]

**33. **N. Ahmed, M. P. Lavery, H. Huang, G. Xie, Y. Ren, Y. Yan, and A. E. Willner, “Experimental demonstration of obstruction-tolerant free-space transmission of two 50-Gbaud QPSK data channels using Bessel beams carrying orbital angular momentum,” in *European Conference on Optical Communication (ECOC)* (IEEE, 2014), pp. 1–3.

**34. **A. Gatto, M. Tacca, P. Martelli, P. Boffi, and M. Martinelli, “Free-space orbital angular momentum division multiplexing with Bessel beams,” J. Opt. **13**, 064018 (2011). [CrossRef]

**35. **J. Arlt and K. Dholakia, “Generation of high-order Bessel beams by use of an axicon,” Opt. Commun. **177**, 297–301 (2000). [CrossRef]

**36. **C. Paterson and R. Smith, “Higher-order Bessel waves produced by axicon-type computer-generated holograms,” Opt. Commun. **124**, 121–130 (1996). [CrossRef]

**37. **S. Qiong-Ge, Z. Ke-Ya, F. Guang-Yu, L. Zheng-Jun, and L. Shu-Tian, “Generalization and propagation of spiraling Bessel beams with a helical axicon,” Chin. Phys. B **21**, 014208 (2012). [CrossRef]

**38. **M. J. Padgett, F. M. Miatto, M. P. Lavery, A. Zeilinger, and R. W. Boyd, “Divergence of an orbital-angular-momentum-carrying beam upon propagation,” New J. Phys. **17**, 023011 (2015). [CrossRef]

**39. **Y. Yan, G. Xie, M. P. Lavery, H. Huang, N. Ahmed, C. Bao, Y. Ren, Y. Cao, L. Li, and Z. Zhao, “High-capacity millimetre-wave communications with orbital angular momentum multiplexing,” Nat. Commun. **5**, 4876 (2014). [CrossRef]

**40. **Brzobohat, “High quality quasi-Bessel beam generated by round-tip axicon,” Opt. Express **16**, 12688–12700 (2008). [CrossRef]

**41. **A. Trichili, T. Mhlanga, Y. Ismail, F. S. Roux, M. McLaren, M. Zghal, and A. Forbes, “Detection of Bessel beams with digital axicons,” Opt. Express **22**, 17553–17560 (2014). [CrossRef]

**42. **A. Dudley, T. Mhlanga, M. Lavery, A. McDonald, F. S. Roux, M. Padgett, and A. Forbes, “Efficient sorting of Bessel beams,” Opt. Express **21**, 165–171 (2013). [CrossRef]

**43. **M. P. Lavery, D. J. Robertson, G. C. Berkhout, G. D. Love, M. J. Padgett, and J. Courtial, “Refractive elements for the measurement of the orbital angular momentum of a single photon,” Opt. Express **20**, 2110–2115 (2012). [CrossRef]

**44. **G. A. Siviloglou and D. N. Christodoulides, “Accelerating finite energy Airy beams,” Opt. Lett. **32**, 979–981 (2007). [CrossRef]

**45. **J. Zhou, Y. Liu, Y. Ke, H. Luo, and S. Wen, “Generation of Airy vortex and Airy vector beams based on the modulation of dynamic and geometric phases,” Opt. Lett. **40**, 3193–3196 (2015). [CrossRef]