
This section describes the applicability of THz PTDH for numerical and experimental studies of achromatic vortex beam modulators. We will demonstrate the capabilities of the PTDH on the example of multilayer quartz wave plates with a certain achromaticity in wavefront retardation in the broad THz range, which have been mentioned in the first paper^{12} of this paired set. However, we note that we propose our technique not only as a method for the quartz platebased beam shapers design, but for any componentbased beam shapers, with any known spectral amplitude and phase characteristics, either full or partially achromatic, such as those reviewed in the first paper of this paired set^{12}. The eligible components can exploit the total internal reflection principle, liquid crystals, metasurfaces and on others principles (such as twisted fibers^{34}, if they are adapted for the THz range). The only necessary condition for the element to be taken into account is the spectral dependence of its action on the incoming beam presented either in the form of spatial amplitude and phase mask, or, for polarising elements, in the form of Jones or Müller matrices.

In the framework of the Jones formalism, linearly polarised along the
$ x $ axis broadband THz field propagating orthogonal to the director axis of a birefringent medium can be represented by the following Jones vector:$$ {\bf{G}}_{{\rm{in}}}(x,y,\nu) = \left( \begin{array}{l} {G_{{\rm{in}}_x}(x,y,\nu)}\\ {0} \end{array} \right) $$ (1) and the action of an optical element onto the electric field vector is given by the Jones matrix.
$$ {J(\nu)} = \left( {\begin{array}{*{20}{c}} A(\nu)&B(\nu)\\ C(\nu)&D(\nu) \end{array}} \right) $$ (2) resulting in a certain frequencydependent change of the broadband THz field:
$$ \left( \begin{array}{l} {G_{{\rm{out}}_x}(x,y,\nu)}\\ {G_{{\rm{out}}_y}(x,y,\nu)} \end{array} \right) = \left( {\begin{array}{*{10}{c}} A(\nu) & B(\nu)\\ C(\nu) & D(\nu) \end{array}} \right) \cdot \left( \begin{array}{l} {G_{{\rm{in}}_x}(x,y,\nu)}\\ {0} \end{array} \right) $$ (3) One of the simplest and most developed approaches to design of achromatic wave plates essential for broadband THz beam manipulation is to stack
$ N $ plates of birefringent material and adjust their thicknesses and orientation to even out the introduced phase retardation of the ordinary wave relative to the extraordinary wave for all spectral components of interest (Fig. 2a). The phase delay introduced by a single$ j $ th birefringent layer, coorientated with the incoming wave polarisation, can be written as:$$ \begin{array}{*{20}{l}} \delta_{j}(\nu,h_j) = {\varphi _x}  {\varphi _y} = 2\pi h_j\left( {{n_{{\rm{e}}}}  {n_o}} \right)\dfrac{{\nu}}{c} \end{array} $$ (4) where
$ n_{\rm{o}} $ and$ n_{\rm{e}} $ are refractive indices for ordinary and extraordinary waves, respectively;$ h_j $ is the thickness of the$ j $ th plate.The influence of each birefringent layer acting as a wave plate in the multilayer achromatic plate on transmitted broadband radiation can be described by the frequencydependent Jones matrix^{35}:
$$ \begin{array}{*{20}{l}} {J_j(\nu)} = \left( {\begin{array}{*{20}{c}} A_j(\nu)&B_j(\nu)\\ B_j^*(\nu)&A_j^*(\nu) \end{array}} \right) \end{array} $$ (5) where
$$ \left\{\begin{array}{l} A_{j}(v)=\cos \left(\dfrac{\delta_{j}(v)}{2}\right)+i \cos \left(2 \theta_{j}\right) \cdot \sin \left(\dfrac{\delta_{j}(v)}{2}\right) \\ B_{j}(v)=i \sin \left(2 \theta_{j}\right) \cdot \sin \left(\dfrac{\delta_{j}(v)}{2}\right) \end{array}\right.$$ (6) Here
$ \theta $ is the angle of the$ j $ th layer optical axis orientation with respect to the$ x $ axis. The resulting matrix of the composite wave plate can be obtained by multiplying the matrices$ J_j $ of each layer in reverse order:$$ J_{a} = \prod\limits_{j=N}^{1}J_j $$ (7) Adjusting the orientation
$ \theta_j $ and thickness$ h_j $ , it is possible to provide any relatively uniform achromatic phase delay$ \delta_a(\nu)\approx{\rm{constant}} $ in a multilayer structure for a relatively broad frequency range^{35}, thus allowing for relatively achromatic wave plate operation of the resulting stack. If the target wave plate is a QWP, the necessary retardation is$ \delta_a ={\pi }/{2} $ in the designed achromaticity interval. Similarly, for the halfwave plate (HWP) the target delay will be$ \delta_a ={\pi} $ .For the designed assembly to work as a qplate, the director of an HWP must be sequentially changed in a circular roundabout of the optical axis. In practise, this is usually achieved by splitting the wave plate into a finite number of sectors. Often, upon manufacture, it is not always possible to reproduce a smooth rotation of the optical axis of a qplate. In such cases, qplate is divided into a discrete number of sectors, each having its own discrete angle
$ \psi $ of the optical axis following the abovementioned rule. This approach allows the introduction of a constant phase delay in each segment$ \Delta \varphi =\pm ({2q\psi }/{M}) $ , where$ M $ is a number of segments selected in such a way that the value of the resulting phase shift throughout the entire plate is equal to$ \varphi = 2\pi \mathsf{L} $ , where$ \mathsf{L} $ is a vortex topological charge. In this case, the dependence$ \theta \left( \psi \right) $ on each mth segment will take the form of a piecewise smooth function, representable by a set of functions in the following form:$$ \begin{array}{*{20}{l}} \Phi \left( \psi \right)=\left\{ \begin{aligned} qm\Delta ,\;\; &\psi \in \left( m\Delta , \left( m+1 \right)\Delta \right) \\ 0, \;\; &\psi \notin \left( m\Delta , \left( m+1 \right)\Delta \right) \\ \end{aligned} \right. \end{array} $$ (8) where
$ \Delta ={2\pi }/{M} $ is the angular measure of each sector,$m=0,1,2,...M1 $ . Thus, the resulting director’s rotation with respect to the plate can be written as the following sum:$$ \begin{array}{*{20}{l}} \theta \left( \psi \right)=\displaystyle\sum\limits_{m=0}^{M1}{\Phi \left( \psi \right)} \end{array} $$ (9) Further extrapolating the idea of designing multilayered wave plates introducing spectrallyisotropic response for THz pulses, we propose the design of an achromatic qplate created from eight sectors of a quartz birefringent achromatic HWP (
$ \delta =\pi $ ) as shown in Fig. 3c in the first paper^{12} of this paired set. Being located in the type I beam shaping scheme (Figs. 1a−d) this plate will introduce uniform topological charge across its whole operational spectrum, and the beam after passing such wave plate becomes a BUTCH beam. Schematic examples of such 8sector qplates for$ \mathsf{L}=1 $ and$ \mathsf{L}=2 $ topological charges and the corresponding cuts of the parent HWPs are shown in Fig. 2b, c.Fig. 3 Frequencyresolved characteristics of the THz field formed by sets of achromatic components arranged in Icx (rows 17) and Ix (rows 12, 89). Rows 13 are the phase differences
$ \Delta \varphi $ between transverse THz field polarisation components, directly after the propagation QWP (row 1), qplate (row 2) and in the detection plane$ z_{\sum} $ for Icxtype (row 3); rows 49 are spectral density modulus (amplitude) (rows 4, 6 and 8) and the phase (rows 5, 7 and 9) images in the detection plane$ z_{\sum} $ for Icx (rows 4, 5 for$ x $  and 6, 7 for$ y $  polarisation) and Ix (rows 89) types. Top insets: the dependencies of the phase shift$ \Delta\varphi $ on the THz radiation frequency introduced by the QWP (left graph) and HWP (right graph). 
Our group designed 7layer achromatic wave plates (Fig. 2) providing a reasonably flat response in the frequency range 0.4 – 1.4 THz. To calculate design, we took the data for refractive indices used in the Ref.^{35}, and minimised the following criteria for QWP and HWP correspondingly:
$$ \begin{aligned}&\displaystyle\sum_{\nu=0.4 {\rm{ THz}}}^{1.4 {\rm{ THz}}}{\\delta(\nu)\frac{\pi}{2}\ + \\Delta\varphi(\nu)\frac{\pi}{2}\ + \\varepsilon(\nu)1\} \\&\displaystyle\sum_{\nu=0.4 {\rm{ THz}}}^{1.4 {\rm{ THz}}}{\\delta(\nu)\pi\ + \\Delta\varphi(\nu)\pi\ + \\varepsilon(\nu)\} \end{aligned} $$ (10) Here
$$ \begin{aligned} &\delta(\nu) = 2\tan^{1}\sqrt{\dfrac{\left{\rm{Im}}\; A\right^2+\left{\rm{Im}}\; B\right^2}{\left{\rm{Re}}\; A\right^2+\left{\rm{Re}}\; B\right^2}}, \; \Delta\varphi=\varphi_y\varphi_x, {\rm{and}} \\ &\varepsilon = \frac{G_x^2 + G_y^2  \sqrt{G_x^4 + G_y^4 + 2G_x^2G_y^2 \cos{2\Delta \varphi}}}{G_x^2 + G_y^2 + \sqrt{G_x^4 + G_y^4 + 2G_x^2G_y^2\cos{2\Delta\varphi}}}. \label{eq:epsilon} \end{aligned} $$ are the chromatic retardation, direct phase delay and the ellipticity, correspondingly.

To track in detail the changes that occur in the spatial frequency structure of the broadband THz vector and/or vortex field, generated by the combination of wave plates and qplates, introduced in Figs. 1a−d, we will use the THz PTDH for a stepwise numerical simulation of the wavefront propagation process^{36}, and supplement it with the effect of the wave plates.
Ideally, all beam shaping components must be stacked together to minimise diffraction and Fresnel losses. However, in practise, it becomes a difficult task. Thus, to avoid superimposition of the Fresnel reflections of the timedomain signal from the wave plateair interfaces, all components of the beam shaper must be separated by a distance
$ z_{\varnothing} > 0.5\cdot c \cdot T_{w} $ , where$ T_{w} $ is the time interval of observation.In PTDH, the equations of scalar diffraction theory are solved for each spectral and polarisation component (see Eqs. (1−3), (6−11) in work^{6}, Eqs. (1−10) in paper^{37}). While detailed description of the method can be found elsewhere^{22,3739}, for the sake of simplicity, the change of the broadband scalar field during propagation to the optical path
$ z $ can be denoted with propagation operator$ {\cal{F}}_z $ . It should be noted that in the propagation process from the source to the first modulator element, the broadband THz field is noticeably transformed^{23}, and therefore the evolution of the field at the distance$ z_{{\rm{in}}} $ must also be taken into account by applying the operator$ {\cal{F}}_{z_{{\rm{in}}}} $ to the field at the source to obtain the vector field in the input plane of the first modulator element:$$ \begin{array}{*{20}{l}} {\bf{G}}_{{\rm{in}}}(x,y,z_{{\rm{in}}},\nu) = {\cal{F}}_{z_{{\rm{in}}}} \left( {\bf{G}}_{{\rm{0}}}(x,y,0,\nu)\right) \end{array} $$ (11) Then, the process of propagation of the THz incoming vector field
$ {\bf{G}}_{{\rm{in}}}(x,y,z_{{\rm{in}}},\nu) $ through the Icxtype beam shaper (shown in Fig. 1a) can be described as:$$ \begin{array}{*{20}{l}} \begin{split} {\bf{G}}_{{\rm{out}}}&(x,y,z_{_{\sum}},\nu) = \\ &= {\cal{F}}_{z_{{\rm{out}}}+z_{{\rm{Icx}}}+z_{\varnothing}}\left(J_{q}{\cal{F}}_{z_{\varnothing}}\left( J_{\frac{\lambda }{4}}{\bf{G}}_{{\rm{in}}}(x,y,z_{{\rm{in}}},\nu)\right)\right) \end{split} \end{array} $$ (12) Here,
$ J_q $ and$ J_{{\lambda }/{4}} $ are the Jones matrices of qplate and QWP, correspondingly;$ z_{_{\sum}} = z_{{\rm{out}}} + z_{{\rm{Icx}}} + 2z_{\varnothing} + z_{{\rm{in}}} $ is the aggregate optical beam path from THz source to the observation plane;$ z_{\varnothing} $ is the distance between wave plates in the modulator;$ z_{{\rm{Icx}}} = h_{{\lambda }/{4}}/n_{\rm{e}} + h_q/n_{\rm{e}} $ is the optical beam path inside achromatic wave plates. For simplicity, we take into account the optical beam path in the wave plates after the matrix transformation of the field, assuming that the change in the structure of the beam when propagating inside these wave plates is insignificant in comparison with the aggregate optical beam path ($ z_{{\rm{Icx}}} \ll z_{\sum} $ ).Similarly, for the type IIx beam shaper, the field at the output can be written as:
$$ \begin{array}{*{20}{l}} \begin{split} {\bf{G}}&_{{\rm{out}}}(x,y,z_{_{\sum}},\nu) =\\ &= {\cal{F}}_{z_{{\rm{out}}}+z_{{\rm{IIx}}}}\left(J_{p}{\cal{F}}_{z_{\varnothing}}\left(J_{\frac{\lambda }{4}}{\cal{F}}_{z_{\varnothing}}\left( J_q{\bf{G}}_{{\rm{in}}}(x,y,z_{{\rm{in}}},\nu)\right)\right)\right) \end{split} \end{array} $$ (13) with
$ z_{{\rm{IIx}}} = h_q/n_{\rm{e}} + h_{{\lambda }/{4}}/n_{\rm{e}} + z_p $ the optical thickness of type IIx modulator components. Here$ z_p \rightarrow 0 $ is the polariser optical thickness and$ J_{p} $ is the Jones matrix of the linear polariser.Thus, by varying the distance
$ z_{{\rm{out}}} $ , one can numerically track the evolution of the formed complexly structured beams during their further propagation in free space. These possibilities of numerical research, after the beam shaper manufacturing, can be supplemented by the potential of physical experimental research. Indeed, having an electrooptical detection module^{40} available or fabricating it using additive manufacturing methods based on open source 3D models^{25,26}, it is possible to record the spatiotemporal dependence of the THz field strength in the diffraction zone. Then, by excluding a possible measurement system response (e.g. excluding noise^{37,41}), one can perform timereversal of broadband THz wavefront and track the reverse dynamics of the pulse propagation, as was repeatedly implemented^{20,23,42}. 
When designing beam modulators, a number of important features should be taken into account. There are 4 factors that should be considered: (i) the limited spectral bandwidth of anisotropic phase delay in broadband achromatic wave plates imposes a restriction on the spectral interval at which the homogeneous topological charge is formed; (ii) the spectrallyresolved ellipticity of broadband polarised radiation, which characterises achromatic QWPs in the composition of the modulator affects the azimuthal homogeneity in the spectrallyresolved amplitude distribution of the formed beam; (iii) qplate introduces the spiral phase shift in the modal distribution of input radiation which leads to the appearance of a singular point, which, in turn, is the cause of an additional diffraction perturbation; (iv) the segmented qplate with eight sectors forms the nonsmooth phase distribution with eight phase steps, the traces of which, during the wavefront propagation, also appear in amplitude.
Design of broadband terahertz vector and vortex beams: II. Holographic assessment
 Light: Advanced Manufacturing 3, Article number: (2022)
 Received: 18 September 2021
 Revised: 13 July 2022
 Accepted: 13 July 2022 Published online: 02 August 2022
doi: https://doi.org/10.37188/lam.2022.044
Abstract:
In this paper, we demonstrate the capabilities of the terahertz pulse timedomain holography in visualisation, simulation, and assessment of broadband THz vortex beam formation dynamics upon its shaping by elements of beam converter, and further propagation and manipulation. By adding Jones matrix formalism to describe broadband optical elements, we highlight the differences in the spatiospectral and spatiotemporal structure of the formed vortex and vector beams dependence on the modulator used and visualise their modal features. The influence of diffraction field structure from each element in the broadband vortex modulator is revealed in numerical simulation and the formed beams are analysed against the simplified LaguerreGaussian beam model.
Research Summary
Design of broadband terahertz vector and vortex beams: II. Holographic assessment
The formation of broadband terahertz vortex and vector beams attracts increasing attention due to the prospects for using such beams in telecommunications, imaging, and sensing of material properties. One of the most effective ways for THz beam shaping is geometric phase approach based on a spinorbit conversion principle. Nikolay V. Petrov from ITMO University and colleagues solved a problem of assessment of broadband terahertz beams properties by application of terahertz pulse timedomain holography for tuning the beam shaping components to ensure maximally possible matching between the characteristics of formed and desired beam. The technique allows estimating the spatiotemporal and spectral characteristics of the THz field inside the modulator and at its output. As examples of beam shaping devices, combinations of multilayer achromatic quartz wave plates as the main functional elements were considered and analyzed.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article′s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article′s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.