- Email: [email protected]

Contents lists available at ScienceDirect

Optics Communications journal homepage: www.elsevier.com/locate/optcom

Theoretical study on 2.52 terahertz beam shaping and polarization conversion based on the transmissive all-dielectric metasurface Zewen Wang ∗, Qi Li, Fei Yan National Key Laboratory of Science and Technology on Tunable Laser, Harbin Institute of Technology, Harbin, Heilongjiang 150081, China

ABSTRACT Based on the physical optics vector diffraction integral formula and the commercial electromagnetic software EastWave, a transmissive all-dielectric square metasurface with a side length of 12.4 mm is designed in this paper. It can simultaneously modulate the phase and polarization of the terahertz wave. The 2.52 terahertz x-polarized normal incident Gaussian beam can be shaped into a radially polarized narrow-width annular Bessel–Gaussian beam (radii of the inner and outer circles are 4.125 mm and 5.5 mm, respectively) at a propagation distance of 23.8 mm away from the metasurface. The metasurface is composed of the poly (4-methyl-1-pentene) substrate and the rectangular silicon resonators with different dimensions and rotation angles (counterclockwise rotation angle of the rectangle long axis as the reference of positive x axis). The relative square error between the amplitude distribution of the shaped beam and the target beam is 7.36%, and the fitting coefficient of them is 96.6%. The diffraction efficiency (defined as the ratio of the total light field energy in the target region to that in the rear plane of the metasurface along the incident beam propagation direction) of the shaped annular region is 75.5%. By focusing the narrow-width annular Bessel–Gaussian beam obtained in this paper, we can get smaller focal spot and longer focal depth, which is of great significance for using terahertz wave in particle acceleration and material processing. In addition, the emergent polarization can be flexibly switched between radially and azimuthally polarized beams by changing the polarization direction of the linearly polarized incident beam.

1. Introduction The cylindrical vector beams have remarkable high-numericalaperture focusing properties [1]. The longitudinal electric field component is strong at the focus of radially polarized beam [2,3], and a tighter focal spot can be obtained than a homogeneously polarized beam [4], which plays an important role in particle acceleration [5] and material processing [6]. Traditional methods for generating radially polarized light include interferometric techniques [7], and the use of a liquid crystal phase modulator [8] or a polarization converter [9]. However, these methods all require additional optical elements that must be aligned, which is difficult to meet the increasing needs for miniaturization and integration. Kitamura et al. [10] show that the depth of focus of the radially polarized Bessel–Gaussian beam can be extended and the spot size made smaller when the intensity distribution of the beam cross section is changed into the narrow-width annular intensity profile. In the latter case the radial component has been strongly suppressed, whereas the longitudinal component dominates the total intensity profile through the focus, thus a smaller focal spot is obtained. The interference of the outermost and innermost rays of the beam axis in the region around the focal point is considered, in the front case, constructive interference take place only at the focal point because the two rays are in phase only at this point. In contrast, the two corresponding rays in the latter case interfere constructively on the beam axis even away from the focal plane, because the incident angles are almost identical, the phases of the rays match well. The latter case

thus possesses a longer depth of focus. They also show that when the width of annular profile is narrow enough, a smaller focal spot size is obtained, compared with the linearly polarized Gaussian beam under the same focusing conditions. The commonly used method to obtain the annular profile is the amplitude modulation of the annular aperture, which will lead to high energy loss. The subwavelength components of the metasurface can flexibly modulate the phase, amplitude and polarization of the electromagnetic wave [11–13]. Compared with traditional optical devices, metasurface devices can be ultra-thin, especially in the terahertz (THz) band [14]. So far, various THz metasurface devices have been reported, such as metalens [14], wave plate [15], beam deflector [16] and beam splitter [17]. Researchers have designed the metasurface which can convert circularly polarized beams into radially polarized beams [18]. However, its efficiency is low due to the inherent ohmic loss of metals. Some multi-layer structures can achieve higher efficiency [19–21], but they need further improvement for massive practical applications due to the problems like alignment of multi-layer patterns in processing [22]. In the THz system, the efficiency can be improved conveniently and effectively by using single-layer transmissive all-dielectric metasurface, which is mainly composed of single-layer dielectric resonator array made of high refractive index materials (such as silicon) and dielectric substrate. It can shape Gaussian beam into rigorous plane wave [23], Bessel–Gaussian beam and vortex beam [24]. However, the transmissive all-dielectric terahertz metasurface that can shape the

∗ Corresponding author. E-mail addresses: [email protected], [email protected] (Z. Wang).

https://doi.org/10.1016/j.optcom.2019.125131 Received 20 September 2019; Received in revised form 17 December 2019; Accepted 17 December 2019 Available online 20 December 2019 0030-4018/© 2019 Elsevier B.V. All rights reserved.

Z. Wang, Q. Li and F. Yan

Optics Communications 460 (2020) 125131

Fig. 1. Schematic diagram of the beam propagation process, amplitude distributions of the incident beam and target beam, as well as the output plane target beam amplitude distribution projections in the polarization direction. (a) Beam propagation process. (b) Incident beam. (c) Target beam. (d) Projection in the x-polarized direction. (e) Projection in the y-polarized direction.

linearly polarized Gaussian beam into the radially polarized narrowwidth annular Bessel–Gaussian beam has not been reported. In this paper, physical optics (PO) vector diffraction integral formula [25] is used to iteratively calculate the continuous phase modulation function corresponding to the metasurface, with the x-polarized ideal Gaussian beam illumination. The waveguide effect of the rectangular silicon resonators with different dimensions can cover 2𝜋 phase modulation and the resonators are arranged on the poly (4-methyl-1-pentene) (TPX) substrate according to the function [26–28]. Then the rotation angle of each resonator is determined to realize the polarization conversion. The metasurface is designed to shape the incident beam into a radially polarized narrow-width annular Bessel–Gaussian beam at a distance of 23.8 mm away from the metasurface, which is of great significance for the applications of the terahertz wave in particle acceleration and material processing. In addition, the emergent polarization can be flexibly switched between radially and azimuthally polarized beams by changing incident light polarization direction [22]. In the second part of this paper, the basic principle of the beam propagation process is introduced. In the third part, a THz transmissive all-dielectric metasurface is designed, and the shaping results are analyzed.

plane center of the metasurface along the incident beam propagation direction. The incident collimated THz wave is shaped into the target beam pattern on the output plane, which is at a distance of 𝑧0 away from the metasurface front plane. A lens with a focal length of f can be placed at the position of output plane to focus the pattern and the focal spot is detected by the detector. The thickness of the metasurface is H. 2.52 THz x-polarized incident ideal Gaussian beam amplitude distribution A(x, y, 0) and the target narrow-width annular Bessel–Gaussian beam amplitude distribution A(x ′ , y ′ , 𝑧0 + H ) on the output plane can be described as: −

𝐴(𝑥, 𝑦, 0) = 𝐶𝑒

𝑥2 +𝑦2 𝜔2 0

−𝛽02

𝐴(𝑥′ , 𝑦′ , 𝑧0 + 𝐻) = 𝐶0 𝑒 { 1 ⋅ 0

√ 𝑥′2 + 𝑦′2 𝐽1 (2𝛽0 ) 𝑟𝑜 √ 𝛼𝑟𝑜 ≤ 𝑥′2 + 𝑦′2 ≤ 𝑟𝑜 else

𝑥′2 +𝑦′2 𝑟2 𝑜

(1)

where C and 𝐶0 are constants. The waist of the incident Gaussian beam is assumed to be located at the rear plane of the metasurface and 𝜔0 is its radius, which is adjusted as half of the metasurface side length 𝐿𝑚 to make the most of the energy being utilized. 𝛽0 is the ratio of target beam outer circle radius 𝑟𝑜 to its waist radius 𝜔1 . 𝐽1 (⋅) denotes the first-order Bessel function of the first kind. 𝛼 is the ratio of the inner radius to the outer circle radius. In this paper, C, 𝐿𝑚 , 𝑧0 , 𝑟𝑜 , 𝛽0 and 𝛼 are chosen as 6.3055, 12.4 mm, 23.8 mm (200 wavelength), 5.5 mm, 1.25

2. Basic principle of the beam propagation process The beam propagation process is shown in Fig. 1(a). The origin of Cartesian coordinate system is assumed to be located at the rear 2

Z. Wang, Q. Li and F. Yan

Optics Communications 460 (2020) 125131

(2) 𝜌′

(x ′2

y ′2 )1∕2

tan 𝜑′

y ′ /x ′ .

where = + and = k is wave vector, 𝛥𝑥 = 𝜌′ cos 𝜑′ − 𝜌 cos 𝜑, 𝛥𝑦 = 𝜌′ sin 𝜑′ − 𝜌 sin 𝜑, R is the distance from the point (𝜌, 𝜑, H ) on the metasurface front plane to the point (𝜌′ , 𝜑′ , 𝑧0 + H ) on the output plane. The total light intensity 𝐼𝑞𝑐 (𝜌′ , 𝜑′ , 𝑧0 +H ) on the output plane can be described as: 2

2

𝐼𝑞𝑐 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻) = |𝐸𝑞 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻)| + |𝐸𝑞𝑧 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻)| ′ (𝜌′ ,

(3)

𝜑′ ,

the radial light field distribution 𝐸𝑞 𝑧0 + H ) on the output plane to be diffracted back to the metasurface front plane can be described as [29]: 𝐸𝑞 ′ (𝜌′ , 𝜑′ , 𝑧0 + 𝐻) = [(1 + 𝜁)𝐴(𝜌′ , 𝜑′ , 𝑧0 + 𝐻) √ ′ ′ − 𝜁 𝐼𝑞𝑐 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻)]𝑒𝑖𝜙𝑞 (𝜌 ,𝜑 ,𝑧0 +𝐻)

(4)

where 𝛷𝑞 (𝜌′ , 𝜑′ , 𝑧0 + H ) is the phase of 𝐸𝑞 (𝜌′ , 𝜑′ , 𝑧0 + H ). 𝜁 is a random constant taken in the range of (0,1). In this paper, 𝜁 is set as 0.5. The radial light field distribution 𝐸𝑞 ′ (𝜌, 𝜑, H ) on the metasurface front plane diffracted backwards by 𝐸𝑞 ′ (𝜌′ , 𝜑′ , 𝑧0 +H ) can be described as: ⎧ ′ (1 − 𝑖𝑘𝑅)𝑧0 𝑖𝑘𝑅 ′ ′ ′ ′ ′ ′ 𝑒 𝜌 𝑑𝜌 𝑑𝜑 ⎪𝐸𝑞 (𝜌, 𝜑, 𝐻) = ∬ 𝐸𝑞 (𝜌 , 𝜑 , 𝑧0 + 𝐻) 2𝜋𝑅3 ⎨ √ ⎪ ′ ⎩𝑅 = (𝛥𝑥′ )2 + (𝛥𝑦′ )2 + 𝑧0 2

Fig. 2. Flow chart of the whole algorithm.

(5)

where 𝛥𝑥′ = 𝜌 cos 𝜑 − 𝜌′ cos 𝜑′ and 𝛥𝑦′ = 𝜌 sin 𝜑 − 𝜌′ sin 𝜑′ . R′ is the distance from the point (𝜌′ , 𝜑′ , 𝑧0 + H ) on the output plane to the point (𝜌, 𝜑, H ) on the metasurface front plane. From 𝐸𝑞 ′ (𝜌, 𝜑, H ) we can get its phase 𝛷𝑞 (𝜌, 𝜑, H ). When it comes to the (q + 1)th iteration, the radial light field distribution 𝐸𝑞 (𝜌, 𝜑, H ) on the metasurface front plane can be described as: 𝐸𝑞 (𝜌, 𝜑, 𝐻) = 𝐴(𝜌, 𝜑, 0)𝑒𝑖𝜙𝑞 (𝜌,𝜑,𝐻)

(6)

finally, the 𝛷𝑞 (𝜌, 𝜑, H ), obtained after a fixed number of iterations, is regarded as the continuous phase modulation function 𝛹 (x, y) corresponding to the metasurface. A transmissive all-dielectric metasurface based on waveguide effect to modulate the phase is designed in this paper. It is composed by the TPX substrate and the subwavelength rectangular silicon resonators with different dimensions and rotation angles 𝜃. Owing to the high refractive index contrast between the silicon cuboid and its surroundings, each cuboid can be considered as a waveguide which is truncated on both sides and operates as a low quality factor Fabry–Perot resonator. Thus, the rectangular cross section of the waveguide leads to different effective refractive indices and phase shifts (𝛷u and 𝛷v), which almost do not vary with rotation angle 𝜃 of the cuboid, of two waveguide modes polarized along the long axis (u axis) and the wide axis (v axis) of the rectangle [12]. The unit cell of square metasurface is composed of TPX substrate and a silicon cuboid with rotation angle 𝜃. When there is a phase difference of 𝜋 between 𝛷u and 𝛷v, each unit cell can behave as a local half-wave plate, the angle of whose optical axis (u axis) and the positive x-axis is 𝜃. Considering a linearly polarized beam normally propagating through an array of half-wave plates (𝜃(x, y) = 𝜑(x, y)/2), the Jones matrix G showing the transformation of the structure is as follows [22]: [ ] cos 𝜑(𝑥, 𝑦) sin 𝜑(𝑥, 𝑦) 𝐺= (7) sin 𝜑(𝑥, 𝑦) − cos 𝜑(𝑥, 𝑦)

Fig. 3. Sketch (left) and top view (right) of the unit cell.

and 0.75, respectively. According to the law of energy conservation, the target beam amplitude distribution A(x ′ , y ′ , 𝑧0 + H ) is determined by the incident beam amplitude distribution A(x, y, 0). Both of them are shown in Fig. 1(b)–(c). When the polarization directions are x and y, the output plane target beam amplitude distribution projections in the polarization direction are shown in Fig. 1(d)–(e), respectively. The continuous phase modulation function 𝛹 (x, y) corresponding to the metasurface can be calculated iteratively by vector diffraction between the front plane of the metasurface and the output plane. The light field initial phase on the front plane of the metasurface is set to be 0. Taking the qth round trip propagation between two planes as an example, we assume that the radial light field is 𝐸𝑞−1 (𝜌, 𝜑, H ) at the point (x, y, H ) on the metasurface front plane, where 𝜌 = (𝑥2 +𝑦2 )1∕2 and tan 𝜑 = y/x. According to PO vector diffraction integral formula [25], the radial light field 𝐸𝑞 (𝜌′ , 𝜑′ , 𝑧0 + H ) and the longitudinal light field 𝐸𝑞𝑧 (𝜌′ , 𝜑′ , 𝑧0 + H ) at the point (x ′ , y ′ , 𝑧0 + H ) on the output plane can be described as:

from Eq. (7) we can see that the incident linearly x-polarized and ypolarized light can be transformed into radially polarized light and azimuthally polarized light after passing through the metasurface, respectively. An arbitrary linearly polarized light can be decomposed into two orthogonal x- and y-polarized lights, so the emergent polarization can in principle be switched between radially and azimuthally polarized beams by changing the polarization direction of the linearly polarized incident beam [22].

⎧ (1 + 𝑖𝑘𝑅)𝑧0 −𝑖𝑘𝑅 ⎪𝐸𝑞 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻) = 𝑒 𝜌𝑑𝜌𝑑𝜑 𝐸 (𝜌, 𝜑, 𝐻) ∬ 𝑞−1 2𝜋𝑅3 ⎪ (1 + 𝑖𝑘𝑅) ⎪𝐸 (𝜌′ , 𝜑′ , 𝑧 + 𝐻) = −𝐸𝑞−1 (𝜌, 𝜑, 𝐻)(cos 𝜑𝛥𝑥 + sin 𝜑𝛥𝑦) 0 ⎨ 𝑞𝑧 ∬ 2𝜋𝑅3 ⎪ −𝑖𝑘𝑅 𝜌𝑑𝜌𝑑𝜑 × 𝑒 ⎪ √ ⎪𝑅 = (𝛥𝑥)2 + (𝛥𝑦)2 + 𝑧0 2 , ⎩ 3

Z. Wang, Q. Li and F. Yan

Optics Communications 460 (2020) 125131

Fig. 4. Phase modulation values and amplitude transmission coefficients of different unit cells. (a) 𝛷𝑢. (b) 𝛷𝑣. (c) 𝛷𝑢 − 𝛷𝑣. (d) Tu. (e) Tv. (f) The relevant data of the chosen unit cells.

The silicon cuboid with dimension D corresponds to phase modulation value 𝛹 (D). The continuous phase modulation function 𝛹 (x, y) corresponding to metasurface is discretized, which means proper dimension 𝐷𝑐 is selected to minimize the value |𝛹 (𝑥, 𝑦) − 𝛹 (𝐷𝑐 )| [26]. In this way, the discrete complex amplitude transmission coefficient of the unit cell at (x, y) can be described as T (𝐷𝑐 )𝑒𝑖𝛹 (𝐷𝑐) . Then the light field distribution 𝐸𝑎 (x, y, H ) of the collimated THz wave which has just passed through the metasurface can be described as: 𝐸𝑎 (𝑥, 𝑦, 𝐻) = 𝐴(𝑥, 𝑦, 0) ⋅ 𝑇 (𝐷𝑐 )𝑒𝑖𝛹 (𝐷𝑐 )

the relative square error (RSE) and fitting coefficient (FC) used to evaluate the shaping effect of the metasurface are denoted as: ∬ [𝐴𝑐 (𝑥′ , 𝑦′ , 𝑧0 + 𝐻) − 𝐴(𝑥′ , 𝑦′ , 𝑧0 + 𝐻)]2 𝑑𝑥′ 𝑑𝑦′ ⎧ ⎪𝑅𝑆𝐸 = ∬ [𝐴(𝑥′ , 𝑦′ , 𝑧0 + 𝐻)]2 𝑑𝑥′ 𝑑𝑦′ ⎪ ⎨ ∬ 𝐴𝑐 (𝑥′ , 𝑦′ , 𝑧0 + 𝐻) ⋅ 𝐴(𝑥′ , 𝑦′ , 𝑧0 + 𝐻)𝑑𝑥′ 𝑑𝑦′ ⎪𝐹 𝐶 = √ ⎪ ∬ [𝐴𝑐 (𝑥′ , 𝑦′ , 𝑧0 + 𝐻)]2 𝑑𝑥′ 𝑑𝑦′ ⋅ ∬ [𝐴(𝑥′ , 𝑦′ , 𝑧0 + 𝐻)]2 𝑑𝑥′ 𝑑𝑦′ ⎩ (11)

(8) E(𝜌′ ,

the diffraction efficiency 𝜂 of the target annular region 𝛺 can be described as:

𝜑′ ,

according to the Eq. (9), the radial light field 𝑧0 + H ) and the longitudinal light field 𝐸𝑧 (𝜌′ , 𝜑′ , 𝑧0 + H ) on the output plane propagated by 𝐸𝑎 (𝜌, 𝜑, H ) can be calculated as: (1 + 𝑖𝑘𝑅)𝑧0 −𝑖𝑘𝑅 ⎧𝐸(𝜌′ , 𝜑′ , 𝑧 + 𝐻) = 𝑒 𝜌𝑑𝜌𝑑𝜑 𝐸 (𝜌, 𝜑, 𝐻) 0 ∬ 𝑎 ⎪ 2𝜋𝑅3 ⎪ ′ ′ ⎨𝐸𝑧 (𝜌 , 𝜑 , 𝑧0 + 𝐻) = ∬ −𝐸𝑎 (𝜌, 𝜑, 𝐻)(cos 𝜑𝛥𝑥 + sin 𝜑𝛥𝑦) ⎪ (1 + 𝑖𝑘𝑅) −𝑖𝑘𝑅 ⎪ × 𝑒 𝜌𝑑𝜌𝑑𝜑 ⎩ 2𝜋𝑅3

𝜂=

∬ 𝐼𝑐 (𝜉, 𝜂, 𝑧0 + 𝐻)𝑑𝜉𝑑𝜂 ∬ [𝐴(𝑥, 𝑦, 0)]2 𝑑𝑥𝑑𝑦

for (𝜉, 𝜂, 𝑧0 + 𝐻) ∈ 𝛺

(12)

where 𝐼𝑐 (𝜉, 𝜂, 𝑧0 +H ) is the total light intensity distribution of 𝛺. The flow chart of the whole algorithm is shown in Fig. 2, and m represents the total number of iterations. The intercoupling between the radial and longitudinal electric field vectors is fully taken into account in this algorithm, which is also applicable to the case of different sampling points and intervals on the input and output planes in the diffraction process.

(9)

the total intensity distribution 𝐼𝑐 (𝜌′ , 𝜑′ , 𝑧0 + H ) and the total amplitude distribution 𝐴𝑐 (𝜌′ , 𝜑′ , 𝑧0 + H ) on the output plane can be expressed as follows: { 𝐼𝑐 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻) = |𝐸(𝜌′ , 𝜑′ , 𝑧0 + 𝐻)|2 + |𝐸𝑧 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻)|2 (10) √ 𝐴𝑐 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻) = 𝐼𝑐 (𝜌′ , 𝜑′ , 𝑧0 + 𝐻)

3. Design of terahertz transmissive all-dielectric metasurface and the analysis of shaping results We use the commercial electromagnetic wave software EastWave [30], which is based on the finite difference time domain (FDTD) 4

Z. Wang, Q. Li and F. Yan

Optics Communications 460 (2020) 125131

Fig. 5. Field plots on the different cross section (x − y cross section: at the position of 30 μm above the TPX–Si interface; y − z cross section: at the position of 𝑥 = 0 μm; x − z cross section: at the position of 𝑦 = 0 μm) of the unit cell 5 with different polarization direction beam illumination. (a) With x-polarized beam illumination. (b) With y-polarized beam illumination.

The thickness H2 of TPX is 118.93 μm. In order to achieve full 2𝜋 phase coverage, according to the relevant principles described in Ref. [33], the period P of the unit cell is chosen to be 55 μm and the height H1 of the rectangular silicon resonator is chosen to be 60 μm. When x-polarized (or y-polarized) light illuminates the unit cell, the variations of its corresponding phase modulation value 𝛷u (or 𝛷v), 𝛷u − 𝛷v (converted to [−180◦ , 180◦ ]) and amplitude transmission coefficient Tu (or Tv) with the lengths Lu and the widths Wv of the

method, to calculate the phase modulation value and amplitude transmission coefficient of the metasurface unit cell (silicon cuboid is located at the center of the unit cell) shown in Fig. 3 (left: sketch, right: top view) at the operating frequency of 2.52 THz, in the case of 𝜃 = 0◦ . The periodic boundary conditions are applied in x- and y-directions, and perfectly matched layers (PMLs) are employed along the propagation direction of the incident light (z-axis) [27]. The refractive index of the silicon is set as 3.418 [31], while that of the TPX substrate is 1.46 [32]. 5

Z. Wang, Q. Li and F. Yan

Optics Communications 460 (2020) 125131

Fig. 6. Variations of the amplitude transmission coefficient and phase modulation value with the imaginary parts of 𝜀TPX and 𝜀Si when y-polarized beam illuminates the unit cell 5. (a) Tv. (b) 𝛷𝑣.

Fig. 7. Simulation results and the proposed metasurface. (a) Continuous phase modulation function. (b) Discrete phase modulation function (8 discrete levels). (c) Unit cells distribution on the metasurface. (d) Diagram of the partial metasurface. Scale bar: 0.275 mm. Table 1 Relevant data of the unit cells 1∼4 corresponding to the Fig. 4.

silicon cuboids are shown in Fig. 4(a)–(e). Therefore, we chose eight different unit cells. The lengths Lu of the silicon cuboids corresponding to the unit cells 1∼4 are 14, 19, 21 and 22 μm, while the widths Wv are 38, 33, 35 and 39 μm, respectively. The unit cells 5∼8 are acquired by rotating the cuboids of the unit cells 1∼4 by an angle of 90◦ counterclockwise as the reference of positive x axis. These unit cells satisfy the conditions: (1) |𝛷𝑢 − 𝛷𝑣| is about 180◦ ; (2) they can realize the phase control of [−180◦ , 180◦ ] for the incident light; (3) the amplitude transmission coefficients remain relatively high at 2.52 THz (the average value is 88.0%). Because the chosen eight unit cells have little difference in amplitude transmission coefficient, all of them can be regarded as local half-wave plates. The data of the eight unit cells corresponding to Fig. 4(a)–(e) are shown in Fig. 4(f), where 𝛷u1 represents the phase modulation value of unit cell 1 and 𝛷u − 𝛷u1 as well as 𝛷u − 𝛷v is converted to [0◦ ,360◦ ]. Relevant data of the unit cells 1∼4 corresponding to Fig. 4 are shown in Table 1. Field plots of the unit cell 5 are shown in Fig. 5, purple and gray patterns mark the outlines of the silicon cuboid and the substrate, respectively. 𝛥z represents the relative distance to the rear plane of metasurface along the incident beam propagation direction. The 2D vector fields on the different cross sections are indicated by the white arrows. From Fig. 5, one will draw the conclusion that the magnetic

Unit cell

1

2

3

4

𝛷𝑢 (degree) 𝛷𝑣 (degree) Tu (%) Tv (%)

−118.8 59.5 95.2 84.2

−72.8 102.3 91.1 89.9

−23.1 147.6 89.0 91.2

26.7 −156.6 90.7 72.8

fields are mainly localized inside the silicon cuboid, which can be seen as the evidence of the waveguide effect [12]. The real parts of TPX permittivity 𝜀TPX and silicon permittivity 𝜀Si are set as 2.1316 and 11.68, respectively. When y-polarized beam illuminates the unit cell 5, the variations of its corresponding amplitude transmission coefficient Tv and phase modulation value 𝛷v with the imaginary parts of 𝜀TPX and 𝜀Si are shown in Fig. 6(a)–(b). With the increase of imaginary parts, Tv gradually decreases. This is because the imaginary part represents the loss term of the material. The larger the imaginary part, the higher the loss. The imaginary part change of 𝜀TPX has greater influence on Tv, because the refractive index of TPX is lower. Compared with Tv, the relative change of 𝛷𝑣 is much smaller and can be neglected. 6

Z. Wang, Q. Li and F. Yan

Optics Communications 460 (2020) 125131

Fig. 8. Output plane amplitude distributions and corresponding projections in the polarization directions. (a) Radial amplitude. (b) Longitudinal amplitude. (c) Total amplitude. (d) Projection in the x polarization direction. (e) Projection in the y polarization direction.

According to the algorithm shown in Fig. 2, the continuous phase modulation function 𝛹 (x, y) obtained by 600 iterations is shown in Fig. 7(a). According to the phase modulation function discretization method mentioned above, Fig. 7(a) is discretized, the result of which is shown in Fig. 7(b). Corresponding to Fig. 7(b), the unit cells distribution on the metasurface is shown in Fig. 7(c). From the color bar of Fig. 7(c), we can see that there are 8 kinds of different unit cells on the metasurface. The diagram of the partial metasurface, which is marked by the black square with an area of 0.66 × 0.66 mm2 at the left bottom corner of Fig. 7(c), is shown in Fig. 7(d). All kinds of unit cells are distributed in that region. According to Eqs. (8) and (9), the calculated radial and longitudinal, as well as the total light field amplitude distributions on the output plane propagated by the collimated THz wave, are shown in Fig. 8(a)– (c), respectively. The output plane amplitude distribution projections in the x and y polarization directions are shown in Fig. 8(d)–(e), respectively. From Fig. 8, it can be seen that the longitudinal light field amplitude is relatively smaller than that of the radial light field, due to the sufficient diffraction distance. The calculated results show that the radial light energy accounts for 99.2% of the total light energy on the output plane. According to Eq. (11), by comparing Fig. 8(c) with Fig. 1(c), we can calculate that RSE = 7.36% and FC = 96.6%. The diffraction efficiency calculated by Eq. (12) is 75.5%. By comparing Fig. 8(d) with Fig. 1(d) and Fig. 8(e) with Fig. 1(e), we calculate that

RSE = 7.24%. Because the unit cells correspond to only eight discrete phase levels and the amplitude transmission coefficients of them are not very close to 1, the simulation results in Fig. 8 are not so ideal compared with the target beam in Fig. 1. The proposed metasurface may be realized in experiments with two steps. The first step is to create the silicon pattern by a sequential process of photolithography and deep reactive ion etching. The second step is to transfer the silicon pattern to the TPX by hot embossing [34–36]. 4. Summary In conclusion, by using the physical optics vector diffraction integral formula and the commercial electromagnetic software, we propose a transmissive all-dielectric square metasurface with a side length of 12.4 mm, which can shape the 2.52 terahertz x-polarized Gaussian beam into a radially polarized narrow-width annular Bessel–Gaussian beam (radii of the inner and outer circles are 4.125 mm and 5.5 mm, respectively) at the propagation distance of 23.8 mm away from the metasurface. The relative square error of the amplitude distributions between the shaped beam and the target beam is 7.36%, while the fitting coefficient is 96.6%. The diffraction efficiency of the shaped annular region is 75.5%. The working frequency of the designed metasurface is 2.52 THz. According to the method in our work, one may design a similar device working at different frequencies, which can 7

Z. Wang, Q. Li and F. Yan

Optics Communications 460 (2020) 125131

simultaneously change the polarization state of the emergent light and shape the incident light into various patterns.

[16] X.Q. Zhang, Z. Tian, W.S. Yue, J.Q. Gu, S. Zhang, J.G. Han, W.L. Zhang, Broadband terahertz wave deflection based on C-shape complex metamaterials with phase discontinuities, Adv. Mater. 25 (33) (2013) 4567–4572. [17] T.M. Niu, W. Withayachumnankul, A. Upadhyay, P. Gutruf, D. Abbott, M. Bhaskaran, S. Sriram, C. Fumeaux, Terahertz reflectarray as a polarizing beam splitter, Opt. Express 22 (13) (2014) 16148–16160. [18] J. Lin, P. Genevet, M.A. Kats, N. Antoniou, F. Capasso, Nanostructured holograms for broadband manipulation of vector beams, Nano Lett. 13 (9) (2013) 4269–4274. [19] C. Pfeiffer, A. Grbic, Cascaded metasurfaces for complete phase and polarization control, Appl. Phys. Lett. 102 (23) (2013) 231116. [20] Z. Wei, Y. Cao, X. Su, Z. Gong, Y. Long, H. Li, Highly efficient beam steering with a transparent metasurface, Opt. Express 21 (9) (2013) 10739–10745. [21] J. Luo, H. Yu, M. Song, Z. Zhang, Highly efficient wavefront manipulation in terahertz based on plasmonic gradient metasurfaces, Opt. Lett. 39 (8) (2014) 2229–2231. [22] F. Zhang, H.L. Yu, J.W. Fang, M. Zhang, S.C. Chen, J. Wang, A.G. He, J.Y. Chen, Efficient generation and tight focusing of radially polarized beam from linearly polarized beam with all-dielectric metasurface, Opt. Express 24 (6) (2016) 6656–6664. [23] Z.W. Wang, Q. Li, Theoretical study on 2.52 terahertz Gaussian beam shaping based on an off-axis transmissive all-dielectric metasurface, OSA Contin. 2 (1) (2019) 201–215. [24] H.F. Zhang, X.Q. Zhang, Q. Xu, Q. Wang, Y.H. Xu, M.G. Wei, Y.F. Li, J.Q. Gu, Z. Tian, C.M. Ouyang, X.X. Zhang, C. Hu, J.G. Han, W.L. Zhang, Polarizationindependent all-silicon dielectric metasurfaces in the terahertz regime, Photonics Res. 6 (1) (2018) 24. [25] T. Hirvonen, J.P.S. Ala-Laurinaho, J. Tuovinen, A.V. Räisänen, A compact antenna test range based on a hologram, IEEE Trans. Antennas Propag. 45 (8) (1997) 1270–1276. [26] M. Khorasaninejad, A.Y. Zhu, C. Roques-Carmes, W.T. Chen, J. Oh, I. Mishra, R.C. Devlin, F. Capasso, Polarization-insensitive metalenses at visible wavelengths, Nano Lett. 16 (11) (2016) 7229–7234. [27] Y.Y. Liang, H.Z. Liu, F.Q. Wang, H.Y. Meng, J.P. Guo, J.F. Li, Z.C. Wei, High-efficiency, near-diffraction limited, dielectric metasurface lenses based on crystalline titanium dioxide at visible wavelengths, Nanomaterials 8 (5) (2018) 288. [28] A. Yariv, Optical Electronics in Modern Communications, Oxford University Press, 1997. [29] Y. Zhao, Y.P. Li, Q.G. Zhou, Vector iterative algorithm for the design of diffractive optical elements applied to uniform illumination, Opt. Lett. 29 (2004) 664–666. [30] The simulation software is EastWave made by Dongjun Information and Technology Co. Ltd. [31] S.H. Ding, Q. Li, R. Yao, Q. Wang, Brewster’s angle method for absorption coefficient measurement of high-resistivity silicon based on CW THz laser, Appl. Phys. B 98 (1) (2010) 119–124. [32] A. Podzorov, G. Gallot, Low-loss polymers for terahertz applications, Appl. Opt. 47 (18) (2008) 3254–3257. [33] J. Zhou, J.J. Wang, K. Guo, F. Shen, Q.F. Zhou, Z.P. Yin, Z.Y. Guo, Highefficiency terahertz polarization devices based on the dielectric metasurface, Superlattices Microstruct. 114 (2018) 75–81. [34] W.S. Lee, K. Kaltenecker, S. Nirantar, W. Withayachumnankul, M. Walther, M. Bhaskaran, B.M. Fischer, S. Sriram, C. Fumeaux, Terahertz near-field imaging of dielectric resonators, Opt. Express 25 (2017) 3756–3764. [35] A. Partanen, J. Väyrynen, S. Hassinen, H. Tuovinen, J. Mutanen, T. Itkonen, P. Silfsten, P. Pääkkönen, M. Kuittinen, K. Mönkkönen, T. Venäläinen, Fabrication of terahertz wire-grid polarizers, Appl. Opt. 51 (2012) 8360–8365. [36] J.J. Wang, K. Guo, Z.Y. Guo, THz filter based on the Si microdisk array, AIP Adv. 9 (4) (2019) 045106.

CRediT authorship contribution statement Zewen Wang: Conceptualization, Methodology, Software, Formal analysis, Investigation, Writing - original draft, Writing - review & editing, Visualization. Qi Li: Conceptualization, Validation, Resources, Writing - review & editing, Supervision, Project administration, Funding acquisition. Fei Yan: Data curation, Visualization. Acknowledgments This work was supported by the National Natural Science Foundation of China (NSFC) under Grant 61377110. The authors thank Yi Zhou of Zhejiang University for his important guidance in metasurface working principle and electromagnetic software simulations. References [1] Q.W. Zhan, Cylindrical vector beams: from mathematical concepts to applications, Adv. Opt. Photonics 1 (1) (2009) 1–57. [2] H.P. Urbach, S.F. Pereira, Field in focus with a maximum longitudinal electric component, Phys. Rev. Lett. 100 (12) (2008) 123904. [3] L.X. Yang, X.S. Xie, S.C. Wang, J.Y. Zhou, Minimized spot of annular radially polarized focusing beam, Opt. Lett. 38 (8) (2013) 1331–1333. [4] R. Dorn, S. Quabis, G. Leuchs, Sharper focus for a radially polarized light beam, Phys. Rev. Lett. 91 (23) (2003) 233901. [5] C. Varin, M. Piché, Acceleration of ultra-relativistic electrons using high-intensity TM01 laser beams, Appl. Phys. B 74 (S1) (2002) S83–S88. [6] O.J. Allegre, W. Perrie, S.P. Edwardson, G. Dearden, K.G. Watkins, Laser microprocessing of steel with radially and azimuthally polarized femtosecond vortex pulses, J. Opt. 14 (8) (2012) 085601. [7] S.C. Tidwell, D.H. Ford, W.D. Kimura, Generating radially polarized beams interferometrically, Appl. Opt. 29 (15) (1990) 2234. [8] G. Miyaji, N. Miyanaga, K. Tsubakimoto, K. Sueda, K. Ohbayashi, Intense longitudinal electric fields generated from transverse electromagnetic waves, Appl. Phys. Lett. 84 (19) (2004) 3855. [9] S. Quabis, R. Dorn, G. Leuchs, Generation of a radially polarized doughnut mode of high quality, Appl. Phys. B 81 (5) (2005) 597–600. [10] K. Kitamura, K. Sakai, S. Noda, Sub-wavelength focal spot with long depth of focus generated by radially polarized, narrow-width annular beam, Opt. Express 18 (5) (2010) 4518–4525. [11] L.X. Liu, X.Q. Zhang, M. Kenney, X.Q. Su, N.N. Xu, C.M. Ouyang, Y.L. Shi, J.G. Han, W.L. Zhang, Broadband metasurfaces with simultaneous control of phase and amplitude, Adv. Mater. 26 (29) (2014) 5031–5036. [12] A. Arbabi, Y. Horie, M. Bagheri, A. Faraon, Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission, Nat. Nanotechnol. 10 (11) (2015) 937–943. [13] Y. Zhou, R. Chen, Y.G. Ma, Design of optical wavelength demultiplexer based on off-axis meta-lens, Opt. Lett. 42 (22) (2017) 4716–4719. [14] X.Y. Jiang, J.S. Ye, J.W. He, X.K. Wang, D. Hu, S.F. Feng, Q. Kan, Y. Zhang, An ultrathin terahertz lens with axial long focal depth based on metasurfaces, Opt. Express 21 (24) (2013) 30030–30038. [15] L.Q. Cong, N.N. Xu, J.Q. Gu, R. Singh, Highly flexible broadband terahertz metamaterial quarter-wave plate, Laser Photonics Rev. 8 (4) (2014) 626–632.

8