Terahertz waveguide-based pulse-shaper system analysis and modeling

In this paper, the transfer function (TF) of a passive waveguide-based terahertz (THz) pulse shaper is calculated using the time domain data provided by the full-wave simulation of the structure. The fractional order of the TF is determined based on the time response resulting from an arbitrary excitation of the proposed pulse shaper. Full-wave electromagnetic numerical analyses are applied to attain the time-domain output data of the helical gold-ribbon dielectric-lined waveguide as the terahertz pulse shaper. To verify the simulation results, the proposed device has been examined using two different numerical methods including the Finite Element Method (FEM) and the Finite Integral Technique (FIT). A good agreement was found between the results of the FIT and the FEM. The use of the system transfer function to analyze the structure is preferable to the full-wave simulation due to execution time savings. Once the TF is determined, one could apply it to subsequent time-domain analyses of the pulse shaper with various inputs.


Introduction
Advances in the field of THz science and technology have made impressive contributions to a variety of research areas including physics, material science (Dorney et al. 2001;Yi et al. 2012;Lopato et al. 2013;Peretti et al. 2018), chemistry (Zhong et al. 2006), and electrical engineering (Hermelo et al. 2017;Grootendorst et al. 2017;Gingras and Cooke 2017; Keren-Zur and Ellenbogen 2019) for several decades. Recently, Veli et al. (2020) presented diffractive layers as a modular pulse shaping network, which is used to shape an arbitrary broadband input terahertz pulse into square pulses with temporal-width using a deep learning technique. After the deep learning-based training of the machine, the thickness profiles of the resulting diffractive layers were determined. Their proposed diffractive networks have been synthesized with various square terahertz pulses with duration from 10.58 ps to 15.69 ps. One of the essential drawbacks of such structures is the bulky size of the pulse shaper (four layers with a spacing of 3 cm ) which is not applicable to incorporation with compact semiconductor terahertz sources. In previous work (Razavizadeh et al. 2020), we demonstrated a tunable THz pulse shaping technique utilizing the positive group velocity dispersion (GVD) feature of a dielectric-lined circular waveguide loaded with a helical graphene ribbon. The results presented in Razavizadeh et al. (2020) indicate feasibility of tunable pulse compression applications by using the graphene ribbon. The current work continues the investigation of a passive fixed pulse-shaper by replacing the helical graphene ribbon with a gold ribbon. Compared with our previous work (Razavizadeh et al. 2020), the pulse shape is greatly improved by reducing the pedestal and ringing. Moreover, in order to validate the numerical analysis, the proposed device has been examined using the FEM and the FIT. The simulated results of the FIT and FEM are in close agreement. There is a need for system identification in engineering (Banos and Gomez 1995), because of its ability to describe phenomena which are not able to be fully explained by analytic approaches in complex problems. The transfer function (TF) modeling is a standard approach like transmission-line modeling and circuit modeling to study the transient behavior of complex electromagnetic problems (Yury et al. 2004(Yury et al. , 2005Zhaona et al. 2021). TF estimation is implemented using a linear model identification method. We iteratively adjust the free model parameters in order to minimize the difference between the full-wave simulation results (such as FIT and FEM) and those of the TF model in response to the input data. Fitting a TF or complex impulse response (CIR) is often used in modern wireless communications in order to characterize wireless channels (Kim et al. 2018;Zhou et al. 2018;O'Hara et al. 2019;Srivastav et al. 2020). The advantage of such a response is to improve the performance of the channel while monitoring the transient responses of the system to various input signals over time-efficient signal processing operations. Examples include the analyses of electromagnetic waves in the presence of complex geometrical boundaries between normal and dispersive materials inducing a transient effect on propagating waves. System identification, which represents the TF, allows a great deal of flexibility in building mathematical models by using measurements of the system's input and output signals in the terahertz spectrum (Zaytsev et al. 2014;Taschin et al. 2017;Huang et al. 2017). Although TFs are generally used in the analysis of systems such as singleinput single-output filters in the fields of signal processing, control theory, and communication theory, it could also be applied to provide comprehensive system models of electromagnetic problems. Figure 1 shows a diagram of the passive waveguide-based terahertz pulse shaper. A gold ribbon of thickness t gold (1 μm), width w (41 μm), and pitch size p (126 μm)is helically wrapped and inserted inside a metalized dielectric tube with inner radius R of 100 μ m and thickness t d of 8.5 μ m. By carefully choosing the geometric parameters and the dielectric material, an optimum condition for an arbitrary dispersion profile, large mode confinement with minimum mismatch, and propagation loss can be achieved. High-density polyethylene (HDPE) is used for the dielectric tube given its good THz transparency (n = 1.524, and = 0.2 cm −1 at 1 THz Zaytsev et al. 2014). The dielectric tube is completely covered by a gold film layer as a metal cladding. Gold is modeled with the Drude parameters as follows: ∞ = 9.1, p = 1.38 × 10 16 rad/s, and = 1.075 × 10 14 Hz (Oral et al. 1983).

Numerical modeling methods for time-domain data acquisition
Numerical techniques are useful in theoretical studies of more challenging electromagnetic problems such as those featuring complex structures where exact analytical solutions are difficult to obtain. In this paper, our proposed device has been analyzed using the FEM and the FIT. Figure 2a, b show the transmittance and absorption coefficients of the proposed device across the whole operating bandwidth (2.14-2.34 THz) for the transverse electric dominant mode (TE 11 ). The phase difference between the input and the output signals of the 2-port device is expressed by the phase of the transmission coefficient ( S 21 ) as depicted in Fig. 2c. The phase of the transmission coefficient, which is numerically calculated using both aforementioned numerical methods, demonstrates the dispersion characteristics of the signal transmission. In fact, it determines how the input signal is delayed when it travels between the input and output ports of the device in the desired frequency spectrum. A passive pulse compressor can be designed in a part of electromagnetic spectra where the normalized group velocity ( V gr ∕c ) profile has a large negative gradient (Piesiewicz et al. 2007;Huang et al. 2017;Macinnes et al. 2006, Bratman et al. 2010, where c denotes the speed of light in free space. V gr is inversely proportional to group delay (GD) when the propagation length (L) is constant as given in Equation (1). On the other hand, GD is defined versus the phase of the transmission coefficient ( ) as stated in Equation (2) (Yang et al. 2008;Ullah et al. 2020).

L(μm) GD(ps)
(2) GD(ps) = − 1 360 (deg) f (THz) As shown in Fig. 3a, the proposed structure shows the normalized group velocity profile with a relatively large negative gradient in the frequency range of f Lo to f Hi which is equal to 2.14-2.34 THz. Furthermore, we provide a comprehensive time-domain analysis that verifies the overall pulse-compression performance of the proposed device. Figure 3b depicts the output envelope waveforms of the proposed device based on a chirpedinput signal obtained using the FIT and the FEM. The chirped-input signal is a y-polarized chirped sine function, as stated in Equation (3): where f Hi is the upper-band frequency of the negative slope of the GV profile (2.34 THz), T 0 is the input pulse width, E 0 is the electric field amplitude, and μ is the chirp factor. The chirp factor for a uniform input pulse spectrum in the negative slope of the GV region from f Lo to f Hi , 2.14-2.34 THz, can be calculated using Equation (4): As the simulation results exhibit, there is good agreement between the results of FIT and FEM over the entire operating bandwidth (2.14-2.34 THz). Theoretically, it is absorption spectra, and c the phase of S 21 of the proposed device using Finite Integral Technique (FIT) and Finite Element Method (FEM) evident that the energy of electromagnetic waves in a dispersive medium just as our proposed structure travels in a GV profile, which is dependent on the frequency as shown in Fig. 3a. Therefore, if such a dispersive structure with negative gradient group velocity profile is excited by using a down-chirped pulse like E y , which is introduced in Equation (3), the frequency spectrum the front pulses with lower GV travel with larger delay of T d,Hi than the tail of the input pulses with delay of T d,Lo as described in Equations (5) and (6) Zhang et al. (2015).
If the input chirped-pulse duration is chosen as (T 0 = T d,Hi − T d,Lo ), one can expect that the tail of the input pulse will overtake the front one. Based on the numerically obtained GV, the optimal traveling length of L 0 is 575μ m, as shown in Fig. 3a. Figure 3a compares the normalized GV and the group delay of the proposed structure with those reported in Razavizadeh et al. (2020). It is clearly seen that the graphene-based pulse shaper in an equal bandwidth around its center frequency, as shown in Fig. 3a, has lower V g than the Au-based one. This effect leads to the longer optimal length (720 μm), and consequently, more transmission loss in the graphene-based pulse shaper. This comparison reveals how the material dispersion contributes to the total chromatic dispersion of the proposed structure. On the other hand, the simplicity is very important in microscale optical the output envelope waveforms of the proposed device based on a chirped-input signal using two different numerical approaches of FIT and FEM, compared with output signal of a zero-biased graphene-based one Razavizadeh et al. (2020) device fabrication. The main challenge in our previous work reported in Razavizadeh et al. (2020) was the design complexity due to the use of graphene and the bias facilities including the SiO 2 spacer and gold-bias electrodes. Therefore, a compact and simple structure able to be incorporated with fibers, signal processors, terahertz time-domain spectroscopy (THz-TDS) systems, and etc. is introduced in this work to relax the structure complexity.
Moreover, Fig. 3b shows that the compressed pulse at the optimal position of the output probe (L 0 ) of the proposed structure has a narrower-shaped pulse with the lower pedestal compared with that of the zero-biased graphene-based structure introduced in Razavizadeh et al. (2020). As shown in Fig. 4, the maximum compressed pulse shape is obtained at the output of a y-polarized electric field, which is located at the distance of 575 μ m from the input reference plane of the proposed device. Figure 4a shows the input pulses just after traveling a length of 50 μ m, about 10% of the optimal length while Fig. 4b-d illustrate the electric field probe outputs located at 525 μ m and 625 μ m, before and after the optimal position (L 0 = 575μm). They clearly show how the probe output signal starts to grow to a peak at the tail position before the optimal position, whereas after passing the optimum length, the temporarily narrowed pulse peak starts to downgrade. The three dimensional (3D) demonstration of the spatial-temporal electric field distribution along the proposed waveguide device can be used to characterize the field confinement performance for future applications such as material detection sensors or time-domain spectroscopy (TDS). Figure 5 gives a comprehensive description of the y-polarized electric field distributions at the optimal probe position in three important moments of the compressed pulse duration, assigned as t A , t B , and t C , which correspond to peak, pedestal, and ringing positions (see Fig. 5a). As depicted in Fig. 5b-d, the field distribution of the peak moment (t A ) ) is significantly more confined than the pedestal (t B ) and the ringing (t C ) moments. This ensures the proposed device can be employed as an efficient pulse shaper for various pulse-based THz applications including fast data communications, TDS, and imaging. Another important conclusion from the electric field distribution would be the mode purity of the proposed structure at the fundamental propagation linear polarized mode, TE 11 , of a circular waveguide. In the remainder of this section, we briefly show the effects of each parameter of the proposed structure on the power transmission and the dispersion performance. Figure 6a shows how the thickness of the Au-ribbon, as helically varying dispersive boundary conditions, influences the GV profile of the proposed structure as well as the insertion loss. The dispersion characteristics of the structure have more dependence on t gold in the lower part of the operating frequency band than elsewhere. Moreover, the return loss of the structure can be tuned effectively within the target operating bandwidth in order to reach the lowest Fig. 4 The output y-polarized E-field probe signals at various positions of the proposed structure a front position (z = 50 μ m), b before, and d after optimal length (575 μm), and finally, at c L 0

Fig. 5
The electric field distributions at the optimum positions and the corresponding peak time positions of b t A , c t B , and d t C as marked on a the probe signal output value with the highest uniformity. Figure 6b shows the internal dielectric coating thickness t d has no significant effect on the GV, especially in the region with negative gradient, while the insertion loss has been kept lower than − 0.4 dB.The dispersion and transmission performance of the proposed structure have significant sensitivity to the most important parameters of the helix which are the gold ribbon width, w, and the pitch size p, which can be seen in Figure 6c-d. Figure 6c illustrates the GV profile is more sensitive to the gold ribbon width below the optimum value, while with a slightly wider ribbon, the GV profile remains unchanged around the center frequency. In addition, for a larger or shorter pitch size, the negative gradient of GV becomes steeper with a down-frequency shift.
In order to ensure that the transient (ringing and pedestal) section has no significant effect on the overall time domain output pulse response, the proposed structure has been simulated in three different fixed lengths of 800 μ m, 1000 μ m, and 1400 μ m, as shown in Fig. 7a-c. The resultant output signals at the optimal position from the input port, L 0 , have been achieved by applying a unique input signal defined with Equation (3). The outcomes show stable responses without any significant deformation which could be used for example for high speed board-toboard interconnecting THz systems (Liu et al. 1996). We can also show that when the structure is excited by an input signal consisting of two sequential chirped pulses, the transient section has no significant effect on the final output compared with that of the single input pulse. This achievement confirms the applicability of the structure to high speed data communications (see Fig. 8).
The next section will describe how to estimate the TF of the proposed structure from the time-domain simulation response.

TF estimation
TF is an important notion in signal processing, which allows one to relate an output signal and an input. This section focuses on the estimation of the general TF of the proposed device. In pulse shaping, knowledge of the TF makes it possible to achieve the temporal output pulse from any arbitrary input pulses. By incorporating the results of the full wave time-domain simulations (presented in the previous section) and the numerical TF estimation approach, it becomes possible to extend the application of the system theory to a wider range of optical problems. The proposed waveguide-based structure is assumed as a linear time invariant (LTI) system which is stable and can be described by a TF as follows: where b 0 , b 1 ,..., b m , a 0 , a 1 ,..., a n are the coefficients to be estimated, and 0 , 1 ,..., m , 1 , 0 ,..., n are the positive real valued exponents. The proposed identification approach is based on the TF estimation using the input/output temporal data which are generated by the fullwave transient simulation of the proposed waveguide under a chirped-input signal. In this work, MATLAB signal processing toolbox was employed to estimate the TF coefficients given in Table 1 based on the obtained time samples of the compressed signal. Figure 9 presents the input-chirped signal and the results of the estimated TF and those obtained by the FEM simulation. It is noteworthy that the coefficients of the numerators and denominators of H(s) with four zeros and nine poles are determined so that the highest similarity . + a n s n Fig. 9 a Chirped-input signal, b reference output signal obtained from the FEM analysis, and c output signal obtained from estimated transfer function model of the proposed device between the FEM output, Fig. 9b, and the reconstructed one by using H(s), Fig. 9c, has been achieved. The accuracy of applying the estimated TF for calculating the output waveform can be verified based on the responses of the proposed structure to different input signals. Figures 10 and 11 show a validation test of the estimated TF H(s) in terms of different input signals ( Fig. 10a Gaussian, Fig. 10b sinusoidal, Fig. 10c chirped-Gaussian, Fig. 10d modulated Rayleigh, Figure 11a-b modulated super Gaussian pulses, and Fig. 11c-d modulated raisedcosine). We briefly introduce all the aforementioned input test waveforms as follows. The chirped Gaussian pulse as a wideband waveform used for validating the estimated TF is stated in the following relation. The first one is a non-chirped sinusoidal pulse while is set to zero and f Hi is replaced with the center operation frequency of f 0 (( f Hi + f Lo )/2) in Equation (3) and the second one is the well-known broadband Gaussian pulse which is defined as Equation (8). where is 2 ps for an appropriate bandwidth compatible with the proposed structure and T 0 is set to 8 ps, the same as previously used in Equation (3).
The chirped Gaussian pulse as a wideband waveform used for validating the estimated TF is stated in the following relation.
where pulse duration ( 0 ), and chirp factor ( μ ) are set to 2.5ps and 15, respectively. When the chirp factor is set to zero, Equation (9) introduces a non-chirped Gaussin pulse, in which the frequency of f 0 should be set at 2.2 THz for better matched bandwidth. All waveforms with time shift of T 0 =5ps are started from time origin. The first order Rayleigh pulse is widely used in ultra-wide band system (Jianxin 2006). To validate the efficacy of the proposed TF, the modulated Rayleigh pulse with center frequency at 2 THz given in Equation (10) is used as the input waveform.
where the pulse parameter ( ) stands for the characteristic time. The small value of corresponds to the narrow waveform in the time domain but to the wide bandwidth in the frequency domain. In our study we set at 1.5ps to have the proper bandwidth. Super Gaussian pulses are commonly used in the modeling of real communications scenarios because Table 1 The coefficients of the numerators and denominators of the estimated transfer function ( m =m, and n = n, where m and n are integer numbers 1, 2, ⋯ ,9) of two reasons: firstly, they usually experience more broadening than the Gaussian pulse and secondly, their shape resembles a real digital pulse of the Non-Return-to-Zero (NRZ) shape (Agrawal 2005). The normalized modulated super Gaussian pulse is described by Equation (11). super-Gaussian with pulse power factor m = 2, and 3, and c, d modulated raised cosine with roll-off factor =1.5, and 0.5, respectively,as well as the corresponding output waveforms obtained from the FEM analysis of the proposed structure where m is an important parameter which has a great influence on the pulse shape. For large values of m, the pulse becomes nearly rectangular with sharp leading and trailing edges. According to Equation (11), the center frequency ( f 0 ), and the pulse duration ( 0 ) are set to 2 THz and 2.5 ps, respectively, in order to have a compatible spectral range with our proposed pulse shaper. A commonly used pulse shape is the raised-cosine pulse which significantly suppresses spectral side lobes and inter carrier interference(ICI). The time domain expression of the modulated raised-cosine pulse is given by Equation (12). A commonly used pulse shape is the raised-cosine pulse, which significantly suppresses spectral side lobes and ICI (Hykins 2014). The time domain expression of the modulated raisedcosine pulse is given by Equation (12).
where is the roll-off factor in the range of 0-1.0. As approaches to zero, the pulse shape becomes more rectangular-like. As shown in Fig. 11c-d, for the roll off factors larger than 1.0, (e.g., =1.5) the raised-cosine pulse seems like a bipolar baseband pulse and a significant difference between the estimated TF output and the FEM results of the proposed device is evident, whereas choosing <1.0 (e.g., = 0.5) or closer to zero (rectangularlike pulse shapes), the estimated TF and the FEM results closely resemble each other. In order to have a compatible spectral range with our proposed pulse shaper, the center frequency ( f 0 ), and the pulse factor (w) are set to 2 THz and, 0.2 THz, respectively. A good agreement between the full-wave results of the FIT simulation and those by applying the equivalent system function, H(s), is observed. As shown in Figs. 10, and 11 this approach has great accuracy and could replace full-wave simulation analysis which is very time consuming. In fact, TF estimation with the features of wideband and memory-efficient analysis could be applied in conjunction with other complicated model estimation techniques such as Transmission-Line (TL) models and narrow band equivalent circuit models.

Conclusions
This paper provides a standard procedure for examining the performance of a waveguidebased terahertz pulse shaper by introducing the TF of the device. The main advantage of this model identification method is that by using a TF, the output of the structure can be easily obtained under different input signals in both time and frequency domains. The TF was calculated based on the input/output time signals of the full-wave simulations of the proposed THz waveguide-based pulse compressor. The zeros and poles of the TF is determined and a close match between the FEM output and the reconstructed output using the estimated TF is obtained. We confirm that all authors contributed to this manuscript. S.M.R. wrote the initial draft of the manuscript as a part of his PhD thesis under advisory of R.S and Z.G.H. and the paper has been reviewed prior to submission by the advisors. S.M.R. provided the concept and simulation design of the study. All authors discussed the results, analyzed the data, and commented on the manuscript. All authors have read and agreed to the published version of the manuscript.