3D-printed THz wave- and phaseplates

: Three-dimensional printing based on fused deposition modeling has been shown to provide a cost-efficient and time-saving tool for fabricating a variety of THz optics for a frequency range of < 0.2 THz. By using a broadband THz source, with a useful spectral range from 0.08 THz to 1.5 THz, we show that 3D-printed waveplates operate well up to 0.6 THz and have bandwidths similar to commercial products. Specifically, we investigate quarter- and half-waveplates, q -plates, and spiral phaseplates. We demonstrate a route to achieve broadband performance, so that 3D-printed waveplates can also be used with broadband, few-cycle THz pulses, for instance, in nonlinear THz spectroscopy or other THz high field applications.


Introduction
Recent progress in the development of THz systems led to a variety of new applications in communications [1][2][3][4], nondestructive testing [5][6][7][8] or particle accelerators [9][10][11][12], to name a few. For THz beam transport and beam manipulation, these systems need different passive THz optical elements, such as mirrors, lenses, beam splitters, waveguides, or waveplates. Beside conventional fabrication methods, such as cutting, milling, polishing, or etching, these elements can also be manufactured using a three-dimensional (3D) printer. Different 3D-printing technologies exist, but the most popular uses fused deposition modeling (FDM). This 3D-printing technique has several advantages: It is time-and cost-efficient, shows good reproducibility, and requires no post-processing. Since the minimum feature size of most 3D-printers is around hundred microns, this fabrication method is best suited for the sub-THz range [13]. Printable materials are usually polymers with a THz refractive index around 1.5, which is favorably low for transmission type optical elements. In addition, a negligible absorption coefficient is essential and is provided, for instance, by polystyrene or by the cyclic olefin copolymer (COC) TOPAS [14,15]. A variety of 3D-printed devices have already been demonstrated, e.g., THz lenses [14,16], gratings [14,16], waveguides [17], axicons [18], phaseplates [19,20], half-waveplates [13] and q-plates [21]. A comprehensive review on different 3D-printed structures for THz applications is found in [15].
Most of these studies were performed in the frequency range between 0.1 THz to 0.2 THz and used continuous wave, narrowband THz sources. To the best of our knowledge, except for Ref. [20] pulsed broadband THz sources have never been used, for instance, to determine the useful bandwidth of 3D-printed THz waveplates. By using such a source, we characterize waveplates in a frequency range from 0.08 THz to 1.5 THz. We consider the most important waveplates for THz applications, namely quarter-waveplates (QWP), half-waveplates (HWP), q-plates and spiral phaseplates (SPP). We show that the operation frequency can be adjusted by a proper scaling of the design and we measure the useful bandwidth of all wave-and phaseplates. Finally, we outline a route toward achromatic waveplates, specifically we demonstrate a QWP design with a twice as large bandwidth compared to the standard QWP. This principle can be extended to other waveplates as well.

Experimental setup
The THz-pulses are generated by a photo-conductive antenna (PCA) gated by a bias voltage and exited by 780 nm, 100 fs laser-pulses at a repetition rate of 80 MHz. We mount a wire grid polarizer to define the THz polarization state to be linear horizontal. The source is imaged with two aspheric lenses with a focal length of 50 mm and 25 mm to the detection plane. In between the two lenses the THz beam has a diameter (1/e 2 ) of about 40 mm and is nearly collimated. Here, we mount the waveplate under investigation. The second lens focuses the modified THz beam and we use a near-field THz scanning microscope to characterize the electric field distribution at the focus as a function of time and transverse position. We cannot measure the electric field directly after the waveplate since the field strength there is too low. However, the measured electric field distribution in the focus is related to the field distribution after the waveplates by a two-dimensional Fourier transformation. Note that we also Fourier transformed the results of numerical simulations and analytical models before comparing them with the measured spatial distributions in the focal plane. The only caveat is that the second lens might modify the THz polarization due to the different Fresnel coefficients for the horizontal and vertical polarization component, but we calculated this effect to be less than 5%.
A detailed description of the near-field scanning unit is given in Supplement 1, Sec. S1, and in Refs. [22,23]. Briefly, the electric field of the THz pulse is measured through electro-optic sampling using ZnTe crystals. The spatial resolution is approximately 20 µm. The detection unit is mounted on a motorized x-, y-and z-stage, such that the electric field can be raster-scanned. By using either a ⟨110⟩-cut or a ⟨100⟩-cut ZnTe crystal we can measure the horizontal, vertical, and longitudinal components of the THz electric field as a function of transverse spatial position and cover a frequency range from 0.08 THz to 1.5 THz. Note that due to imperfections of the coating and of the ZnTe-crystals, the detection efficiency varies for the different crystal orientations, and we estimate these variation to be smaller than 10 %. Based on the statistical analysis of multiple measurements we estimate the error in the measured phase to be on the order of 0.1 rad in the frequency range of interest. In the following we consider these amplitude and phase errors in the error propagation.

Wave-and phaseplate design and fabrication
In order to manufacture a birefringent metamaterial, out of which the waveplates were fabricated, we consider a one-dimensional array of alternating layers with a refractive index of n 1 and n 2 with a corresponding layer thickness of L 1 and L 2 . The effective refractive indices n e and n o for an electric field polarization parallel and perpendicular to the layers are given by [24,25] where Λ = L 1 + L 2 is the period of the structure, λ is the free-space wavelength, and η 1 = L 1 /Λ and η 2 = L 2 /Λ are the filling factors of the individual components. We defined n 2 ∥ := η 1 n 2 1 + η 2 n 2 2 and n 2 (1) is only strictly valid for λ>n 2 Λ [26] assuming n 2 >n 1 . In the limit of λ ≫ Λ, the indices of refraction asymptotically approach n e ≈ n ∥ and n o ≈ n ⊥ . For lower wavelengths, diffraction from the periodic one-dimensional grating structure degrades the performance of the waveplate.
Prior to fabrication, all waveplate designs were optimized by time-domain simulations performed with CST [27]. The simulation volume contained a single unit cell of the grating and we used periodic boundary conditions in transverse directions. The incident THz pulse was modeled as a plane wave and open boundary conditions were used in the propagation direction. With the help of the simulations the design parameters of all waveplates, i.e., grating period, height and filling factor, were optimized within the limitations of the printing process. Specifically, we fixed the wall thickness to a value determined by the nozzle of our 3D-printer. Then, we ran parametric sweeps for the period and the structure height and selected the best combination.
The wave-and phaseplates were fabricated using standard FDM technology with an Ultimaker 2+ 3D-printer. We use commercially available COC-filaments, which have a measured refractive index of 1.53 and an absorption coefficient of less than 1 cm −1 for frequencies between 0.2 THz to 0.8 THz [14]. The waveplates are based on a grating-like alternating air (n 1 = 1) and filament (n 2 = 1.53) structure with the design parameters calculated from Eq. (1). All waveplates were printed in less than 12 hours using a nozzle diameter of 0.25 mm and a layer height of 0.1 mm. The nozzle temperature was set to 210 • C and the build plate temperature to 70 • C. All the waveplates have a 3 mm thick 3D-printed cylindrical support plate with a diameter of 50 mm. In order to minimize imperfections from 3D-printing we selected the wall thickness for all waveplates such that walls are printed as single lines from the extruder. Measurements revealed a wall thickness of (0.22 ± 0.02) mm, which is slightly smaller than the nozzle diameter. The error of 20 µm is due to thickness variations in the printing process. In order to minimize their influence we used, whenever possible, a filling factor well below 0.5. As a result the 3D-printed structure shows a birefringence of approximately 0.12, which is more than two times higher than that of standard THz waveplates made out of quartz [28]. This is of special interest for the low THz frequency range, where standard waveplates tend to become rather bulky. From the frequency dependent birefringence ∆n(f ) we calculate the structure height to h = δ  First, we printed QWPs and HWPs with a design similar to the one proposed in Ref. [13], but scaled to the more than two times higher design frequency of 0.35 THz. The grating structure has a period of 0.7 mm and the grating height is 1.9 mm for a QWP and 3.8 mm for a HWP. To evaluate the influence of the grating period, we manufactured an additional HWP with a 1.1 mm period while keeping the grating height constant. Moreover, to demonstrate the tunability we fabricated a QWP and a HWP both with a higher design frequency of 0.5 THz. Here, the grating period was reduced to 0.4 mm and the grating height to 1.5 mm (QWP) and 3.0 mm (HWP). Next, we printed q-plates, which locally act as HWPs whereby the orientation of the fast axis changes as a function of azimuthal angle. For a q-plate with a topological charge of q = 0.5, the angle of the fast axis changes as α(φ) = 0.5φ + α 0 , where φ is the azimuthal angle and α 0 is an offset corresponding to the plate orientation. The layout of the grating structure consists of a set of profile curves in the xy-plane each given by y(x) = ± √︁ y 0 (y 0 + 2x), where y 0 is the intercept of the considered curve with the y-axis. Note that the grating period and the filling factor change depending on the distance to the origin. The profile curves are chosen such that the grating period is 1.1 mm at a radius of 0.8 cm. Closer to the origin, the grating period is smaller, whereas further away it is larger. The grating height is fixed to 3.8 mm.
Spiral phaseplates, which change the orbital angular momentum of the THz beam, do not require a birefringent material, but a phase retardance that decreases or increases linearly with azimuthal angle. For a SPP with a topological charge of ±1 at the design frequency f , the step height h at an azimuthal angle of zero is given by [29,30] where ∆n(f ) is the frequency dependent refractive index difference between the plate material and its surroundings. As an example, we consider a topological charge of ±1 at 0.3 THz, which corresponds to a step height of 1.9 mm. As demonstrated in [18] also higher topological charges can be achieved by using integer multiples of the step height h or multiple equidistant steps along the azimuthal angle. Ultimately, the step height is limited by the minimum 3D-printed layer thickness, which is typically as low as 60 µm.

Results
In this section, we analyze the performance of a QWP, HWP, q-plate, and SPP as a function of frequency, compare the results to an analytical model as well as numerical simulations, and extract the useful spectral bandwidth. Even though the setup allows to record data up to 1.5 THz, we show data only up to 0.6 THz, because the onset of diffraction above the critical frequency f c = c/Λ leads to meaningless results. In addition, diffraction losses are not properly captured by the CST simulations, because we use periodic boundary conditions. The analytical model is based on 4f -imaging of a linearly polarized Maxwell-Gaussian beam [31,32], where the polarization or the phase in the Fourier plane (2f ) is modulated in accordance with an ideal wave-or phaseplate performance. A design frequency around 0.35 THz is a good compromise between the 3D-printer resolution and the spatial resolution of the THz near-field scanning setup. However, for one example, i.e. a QWP, we demonstrate that the design frequency can be adjusted by a scaling of the design parameters. Finally, we consider an achromatic QWP design based on a two layer structure.

Quarter-waveplate
Generally, a QWP produces an elliptical state of polarization with a left-or right-handed symmetry. To verify the proper functioning of the 3D-printed QWP, we measured the horizontal and the vertical polarization component for different waveplate rotation angles. In the following, a rotation angle of 0 • corresponds to the fast axis pointing in (1, 0, 0)-direction. Figure 2(a) shows the birefringence as function of frequency. We compare the experimental results with CST-simulations as well as the analytic model based on Eq. (1). While the agreement between experiment and CST simulations is close to perfect for the entire frequency range, the analytical model starts deviating above 0.45 THz, where it is no longer valid since it neglects diffraction at the line structures. A key-parameter to characterize the performance of the QWP is the phase retardance between the fast and the slow axis. Figure 2(b) shows the retardance normalized to π as function of frequency with the blue circles referring to the experimental results and the dashed curve corresponding to the analytic model assuming a standard QWP with a frequency independent birefringence. The grey region indicates the frequency range in which the retardance deviates less than ±10% from the ideal value of π/2. For the experimental results we find a bandwidth of approximately 60 GHz (i.e., a relative bandwidth of 17%), while the analytical model for a standard QWP predicts 70 GHz (i.e, a relative bandwidth of 20%). Thus, the bandwidth of the 3D-printed QWPs is similar to that of a standard commercial QWP. At the design frequency we found that approximately 80% of the incident power is transmitted and we attribute 9% of the losses to Fresnel reflections at the two interfaces and 11% to material absorption [14]. For comparison, the transmitted power for a standard waveplate made of quartz with n = 2.1 [28] is 76%. To study the frequency dependent waveplate performance it is useful to consider the Stokes parameters. Note that we can calculate the Stokes parameters from only two measurements, because we have access to the x-and y-component of the electric field. The normalized Stokes parameters are shown in Fig. 3(a)-(c) where the solid black curves correspond to simulation results. We find excellent agreement between experiment and simulations. At the design frequency of 0.35 THz the normalized Stokes parameters of an ideal QWP are S 1 /S 0 = 0, S 2 /S 0 = 0, and S 3 /S 0 = ±1. For a QWP rotation angle of 45 • the polarization is left-handed circular (using the source-convention) and we find S 1 /S 0 = −0.23 ± 0.15, S 2 /S 0 = 0.20 ± 0.13 and S 3 /S 0 = −0.95 ± 0.04. For a 135 • rotation angle the polarization is right-handed circular with S 1 /S 0 = −0.10 ± 0.16, S 2 /S 0 = −0.19 ± 0.14 and S 3 /S 0 = 0.98 ± 0.03. To change the design frequency, we simply have to scale the design parameters, which we demonstrate by fabrication of a QWP for 0.5 THz. When we combine the waveplate with a bandpass filter also centered at 0.5 THz and with similar bandwidth, we expect perfect circularly polarized THz pulses. At 0.5 THz the filter transmits 80 % of the power and the bandwidth at 1/e 2 is 0.15 THz. Figure 4 shows the measured horizontal and vertical component of the THz electric field for two different orientations of the QWP, namely at 45 • (a) and at 135 • (b). For better visualization, we applied a moving averaging with a 300 fs time window to the measured signals. While the horizontal components are identical for both orientations, the vertical components show a π phase-shift. As a result, the electric field vector as a function of time follows a left-handed or right-handed spiral, as indicated by the red curves. Hence, we can unambiguously distinguish between left (a) and right (b) circular polarization.

Half-waveplate
Next, we consider a HWP and Fig. 5 shows the relative power of the x-(a) and y-polarization (b) as a function of frequency and HWP rotation angle. At rotation angles of 45 • and 135 • and for frequencies around 0.35 THz the horizontally polarized component (x) is minimal, whereas the vertical component (y) is maximal. At lower frequencies, the polarization remains predominantly horizontal, and at higher frequencies, diffraction losses at the HWP structure set in. These results indicate the correct functioning of the HWP with a phase retardance of π at the design frequency. Figure 3(d)-(f) show the normalized Stokes parameters as function of frequency for a HWP rotated to 45 • and 135 • . For comparison we show the simulated results for a rotation angle of 45 • and we find excellent agreement between experimental results and simulations. Around the design frequency of 0.35 THz the polarization is rotated by 90 • and the experimentally determined normalized Stokes parameters are S 1 /S 0 = −0.999 ± 0.001, S 2 /S 0 = 0.016 ± 0.006, and S 3 /S 0 = 0.035 ± 0.006. Note that a HWP performs as a QWP around one half of its design frequency. Indeed, Fig. 3(d)-(f) show normalized Stokes parameters indicative of a QWP around 0.2 THz. This frequency is larger than half of the design frequency, which is explained by the fact that the birefringence increases with frequencies, as shown in Fig. 2(a). The influence of the grating period on the performance of the HWP is of special importance for the design of the q-plates, since they consists of a grating structure with a spatially varying grating period. Therefore, we increased the grating period of the HWP from 0.7 mm to 1.1 mm while keeping all other parameters fixed. For the larger period of 1.1 mm, diffraction losses start already at a frequency around 0.27 THz. Nevertheless, we find a normalized Stokes parameter S 1 /S 0 = −0.987 ± 0.004 at 0.3 THz, which indicates that the HWP performs as expected even if diffraction losses decrease the overall transmission.

q-plate
q-plates, depending on their rotation angle, convert an incident linear polarization into either a radial or an azimuthal polarization. If the symmetry axis of the waveplate is parallel (perpendicular) to the incident linear THz polarization, the beam is converted to a radially (azimuthally) polarized beam. To verify the proper functioning of the q-plate we measure, spatially resolved, all three electric field components. We scan an area of 3 mm-by-3 mm in the xy-plane in steps of 0.1 mm (i.e., 31-by-31 points). Figure 6 shows the field distribution of the radial (a) and the azimuthal (b) polarization at a frequency of 0.325 THz. The color-coded amplitude of the transverse electric field is shown together with a vector field indicating the instantaneous direction. The amplitude shows the typical doughnut-shape spatial profile with a singularity at the origin. The three insets show the amplitude distribution of the horizontal, vertical, and longitudinal electric field components in the transverse plane. The horizontal and vertical field components exhibit the characteristic spatial distribution with two lobes, which are π-phase shifted. The asymmetry of the azimuthally polarized beam profile is most likely due to an unwanted tilt of the q-plate. In Supplement 1, Sec. S4B, we present an analytic model assuming a perfect q-plate in a 4f-imaging system and we find good agreement of the resulting spatial polarization pattern with the experimental results. Moreover, for the radially polarized beam, we observe a strong longitudinal electric field component close to the beam center. A detailed description of the longitudinal field for a radial THz pulse can be found in [33]. For an azimuthal polarization, the longitudinal electric field is negligible, instead there is a strong longitudinal magnetic field component (see Supplement 1, Sec. S4B). To characterize the polarization state as a function of frequency, we project the measured transverse electric vector field distribution E(x i , y j , f ) onto an ideal mode of the same polarization state where the sum runs over all pixels and the amplitude of E 0 is chosen such that a perfect radial or azimuthal polarized mode yields γ(f ) = 1. The global phase of the electric field distribution E(x i , y j , f ) is adjusted for a maximum γ(f ). The beam waist w 0 (f ) is obtained from a fit to the measured amplitude distribution and we found w 0 (f ) = 600 µm 0.325 THz f at the frequency f . Figure 6(c) shows the projected polarization, γ, as function of frequency. A similar performance is found for both q-plate orientations, with the maximum projected polarization of 0.95 at 0.325 THz being close to the ideal value of 1. For frequencies between 0.275 THz and 0.375 THz we find a projected polarization larger than 0.9, which indicates that the q-plate works as intended in this frequency range. At lower and higher frequencies the projected polarization decreases rapidly.

Spiral phaseplate
SPPs are used to generate beams with angular momentum without affecting the state of polarization. Hence, the transmitted THz pulse remains horizontally polarized, but depending on the spatial position the wavefront is retarded. We measured the horizontal electric field component in the transverse plane in an area of 3 mm-by-3 mm in steps of 0.1 mm (i.e., 31-by-31 points). Figure 7 shows the amplitude of the horizontal field component at 0.3 THz (a) and the corresponding phase (b). The phase is shown only when the amplitude at the same pixel is larger that 30% of the maximum amplitude. The amplitude distribution shows a doughnut-shape spatial profile due to a singularity at the beam center. The color-coded phase distribution indicates that the phase retardance decreases linearly by 2π as function of the azimuthal angle around the beam center and exhibits a phase step ∆ϕ at 0 • . Due to the layer resolution of 0.1 mm, which is ten times smaller than the considered wavelength, no discrete phase steps are observed. Figure 7(c) shows the phase step ∆ϕ as function of frequency. The phase was averaged over the radial direction around the amplitude maximum, i.e., between r = 0.3 mm and 0.6 mm, as indicated by the two white circles in Fig. 7(b). The phase step ∆ϕ normalized to 2π corresponds to the topological charge. Note that ∆ϕ has to be an integer multiple of 2π, because the unwrapped phase of electromagnetic fields must be continuous. For comparison, we used the same recipe to calculate the phase step from an analytical model. In both cases, we find discrete steps when plotting the phase step as function of frequency. When we use alternative recipes to extract the phase step, the overall picture remains the same, except the exact frequencies at which the topological charge changes are found to shift to somewhat lower or higher frequencies. On average, we find the first two steps at 0.15 THz and at 0.45 THz. These frequencies corresponds to 0.5 times and 1.5 times the design frequency, respectively. We also tested two different SPPs with a topological charge of ±2 at 0.3 THz and found similar performance.

Achromatic waveplates
For broadband THz applications, for instance using a femtosecond laser pumped LiNbO 3 source, the useful bandwidth of the waveplates presented above is not sufficiently broad. Hence, we aim at increasing the bandwidth by using an achromatic waveplate design. As an example, we discuss a QWP with an achromatic design consisting of two layers where the fast axis of the first layer is aligned with the slow axis of the second layer. Naturally, this design is motivated by commercially available achromatic waveplates. To optimize the performance, we tune the two grating heights as well as the two periods. Hence, we have a four dimensional parameter space, which allows us to optimize the structure for broadband performance. For optimization we use a genetic optimization algorithm (in Matlab) and minimize the cost function ∫ df |φ(f ) − π/2| for f ∈ [0.23 THz, 0.4 THz], where φ(f ) is the phase retardance based on Eq. (1). We impose boundaries to the optimization: 1) The grating period is limited to the range between 300 µm and 675 µm, where the lower limit is due to the 3D-printer resolution and the upper limit is calculated from the onset of higher order diffraction, and to integer multiples of the accuracy of the 3D-printer, i.e., 12.5 µm.
2) The grating height is restricted to values lower than 10 mm, in order to be able to mount the achromatic QWP in a standard rotation mount, and to integer multiples of the minimum layer thickness, i.e., 100 µm. After successful optimization we find the following design parameters: The first grating period is 412.5 µm and has a height of 5.8 mm. The second grating has a grating period of 662.5 µm and a height of 3.5 mm. Between the two layers we include a 1 mm thick spacer layer for mechanical support of the second layer. Note that the achromatic structure was printed as a single part. Figure 8 compares the retardance of the single layer QWP with the optimized achromatic design. The grey area indicates a 10% deviation from the ideal retardance of π/2. The experimental data are in good agreement with the double layer simulation results. The bandwidth of the achromatic design is more than two times larger than that of the single layer design. Diffraction losses limit the useful frequency range up to approximately 0.5 THz. We would like to stress that the bandwidth can be even further increased by using more than two layers and that the same design principles can be applied to HWPs and q-plates. Fig. 8. Measured retardance as function of frequency for the original one layer QWP (blue circles) and for the optimized achromatic design (red stars). The blue circles corresponds to the results already shown in Fig. 2(b) and the green curve corresponds to the simulation of the achromatic design. The grey area indicates a 10 % deviation from the ideal retardance of π/2.

Conclusions
We experimentally evaluated 3D-printed QWP, HWP, q-plate, and SPP for broadband THz applications and compared the results to finite element simulations. We demonstrated that 3Dprinted wave-and phaseplates, using fused deposition modeling, are not limited to a frequency range <0.2 THz, as indicated by previous publications. By using a broadband THz source with a useful spectral range from 0.08 THz to 1.5 THz we showed that 3D-printed waveplates perform as expected up to about 0.6 THz. We found that all waveplates have bandwidths similar to commercial products and we demonstrated a route to achieve broadband performance, so that 3D-printed waveplates can also be used in combination with few-cycle THz pulses. Our results demonstrate that 3D-printing is a cost-efficient and time-saving fabrication method for a variety of narrow-and broadband wave-and phaseplates.