Novel technique for medium permittivity profile reconstruction using THz pulsed spectroscopy

A new method for permittivity profile reconstruction based on terahertz time-domain spectroscopy signal processing is presented. Reconstruction is accomplished in two steps. First, the sample pulse function is reconstructed using sample time-domain reflection data. Low and high frequency noise filtering and the interpolation of the pulse function at low frequencies are then applied. Second, an invariant embedding technique is used to calculate the dielectric permittivity profile based on the sample pulse function. Method was tested experimentally on samples with known permittivity profiles. The algorithm is stable to additive Gaussian white noise as shown using mathematical modeling based on the finite-difference time-domain technique (FDTD). Possible applications of this permittivity profile reconstruction technique are discussed.


Introduction and background
Terahertz (THz) radiation is located between the infrared and the microwave regions of the electromagnetic spectrum, from about to . THz waves have some specific properties, such as high penetration depth into dielectric materials, and the ability to interact with both vibrational and rotational molecular levels [1]. THz radiation is non-ionizing in nature. Because of these properties, there are many potential applications of THz technology [2]- [5].
THz time-domain spectroscopy (TDS) utilizes short pulses of THz radiation with a wideband spectrum (from to ) to measure the THz material parameters of various media. The electric field of a THz pulse is detected with a very high time resolution (about ) after the transmission of the pulse through the sample or the reflection of the pulse from the surface of the sample. By applying a fast Fourier transform to the detected time-domain data, it is possible to analyze the complex transmission or reflection coefficients of the sample in a wide frequency range.
Time-domain data of TDS system can also be used to reconstruct the permittivity profile in a sample, i.e., the dependence of the sample dielectric permittivity on depth. This reconstruction method aids in the study of the internal structure of a sample, known as THz tomography or T-Ray tomography [6]. THz tomography is expected to be useful in medical diagnosis, nondestructive evaluation of construction materials, inspection of art objects, and other applications. This paper describes a new way to solve this inverse problem. The results of implementing the new algorithm and its stability with respect to noise are shown.
Several assumptions concerning media properties are made in the present work. The medium of interest is assumed to be non-dispersive and non-absorbing. Sample optical properties are assumed to be constant in the lateral direction and vary only with depth. The media of interest must possess no optical non-linearity. The medium is assumed to be thick such that multiple reflections in sample layers do not impact the reflected signal data (i.e. the number of pulses in the sample response). Notice, while studying materials with rather high dispersion of optical properties we are reconstructing sample permittivity profile, which corresponds to the average of the sample permittivity in the entire bandwidth of the THz pulse spectrum.
We consider all stages of permittivity profile reconstruction using this algorithm. We discuss all methods for solving the problems of invariant embedding technique generalization to THz spectroscopy.

Algorithm description
Sample dielectric permittivity profile reconstruction is accomplished in two steps: reconstruction of the sample pulse function based on TDS system signals and the reconstruction of the sample permittivity profile based on a sample pulse function . Two signals need to be obtained for sample permittivity profile reconstruction: is the signal reflected from the sample of interest, and is the reference signal reflected from the reference surface with a very high and homogeneous reflectivity in a wide spectral range. A planar gold mirror was used to obtain the reference signal.
The sample pulse function is the medium response to the influence of an incident electromagnetic wave with amplitude time dependence, where is a Dirac delta-function. The sample pulse function is also the result of an inverse Fourier transform applied to the sample amplitude reflectivity ̃ . A technique utilizing a Wiener filter ̃ was applied in previous research [7] to process TDS system data. The formulation of a derived Wiener filter with several modifications is used to solve the deconvolution problem in Fourier-domain and to obtain sample spectral reflectivity ̃ : where is a noise power spectrum model, and is a signal power spectrum model. White Gaussian noise power spectrum as a basis for the function . The function could be easily defined with Gaussian monopulse model of TDS signal. The equation (1) is negligible in a frequency range where the noise power spectrum model exceeds the signal power spectrum model . Conversely, if exceeds , the function ̃ is close to unity. Smooth filter edges are also provided. This filtering procedure allows us to obtain as much data on the sample reflectivity as possible, and contains smoothed filter edges to prohibit the emergence of the Gibbs noise in the sample pulse function [8], [9]. Obtained sample reflectivity depends only on the dielectric permittivity ε(z) and the effective conductivity σ(z) of the sample, and does not depend on the properties of TDS system signals.
The second problem in the pulse function reconstruction procedure is a reflectivity interpolation at low frequencies. To reconstruct a sample permittivity profile with the invariant embedding technique, it is essential to know the sample pulse function spectral components at low frequencies, excluding the zero point of the discrete Fourier-domain. It is possible to interpolate the sample reflectivity ̃ based on known frequency domain information within the ranges [ ] and [ ], and to obtain low frequency data from the interpolation results. It is convenient to use trigonometric interpolation for this purpose [8], [9]. We apply this type of interpolation for the modulus | ̃ | and the phase { ̃ } of the sample complex reflectivity independently. Trigonometric interpolation of the real { ̃ } and imaginary { ̃ } parts of the sample reflectivity could provide good reconstruction results too. After interpolation the sample pulse function reconstruction take place according to the equation: ( ̃ ) where is the inverse Fourier transform operator, is the pulse function containing lowfrequency components.
To determine the dielectric permittivity profile ε(z) of the sample based on the sample pulse function we need to solve the following initial-boundary problem with the sample pulse function as an initial condition: ∫ √ is a time which takes the wavefront to travel through the sample with thickness , ∫ is a normalized spatial coordinate, -is a normalized time coordinate, is a normalized integral kernel corresponding to the sample depth . Note, that is the initial condition for computations. The function is a result of the initial-boundary condition (3) numerical solution. Permittivity profile can be easily determined based on utilizing the following transformation [10]: where is the speed of light in region of space , is the permittivity corresponding to the region of space . The algorithm can be modified to take media absorbance into account. One can use some a priori information about the sample conductivity profile or some approximation of the conductivity distribution during reconstruction. In addition, it is possible to connect the sample conductivity with the sample permittivity utilizing a mathematical model of the complex permittivity ̃ and Kramers-Kronig relation. In this case, the reconstruction of some ̃ parameter profile will take place.

Experimental results
In order to validate this algorithm, several media with known permittivity profiles were studied. These test media were built from a set of flats with known thickness and THz optical properties. Figure 1 shows an example of the reference and the sample signals detected during the study of the test sample build of thick quartz window. The results for an algorithm implementation and the priori data

Algorithm stability in the presence of white noise
The stability of the permittivity profile reconstruction algorithm was considered by mathematical modelling of the algorithm implementation in the presence of additive white Gaussian noise in TDS system signals. Models of layer medium permittivity profiles as well as models of permittivity profiles with slightly smoothed edges were considered. The interaction of terahertz light with matter was examined utilizing FDTD modelling and reflected waveforms were obtained. White Gaussian noise with various standard deviations was added to the waveforms and the permittivity reconstruction procedure was applied. Results of algorithm implementation were compared with the initial model. Dependence of the standard deviation for the profiles on the waveform noise standard deviation for the same function is shown in figure 3. This figure also contains the border between stable 1) and unstable 2) regions. The dependence of the standard deviation of permittivity profile reconstruction on the noise standard deviation is less than 15.0% for rather high SNR of time-domain data.

Conclusions
In the present work, the ability of THz time-domain spectroscopy applied to media dielectric permittivity profile reconstruction was shown. Novel algorithm for solving the inverse scattering problem based on an invariant embedding technique was developed and implemented. The present algorithm can be used in various fields of THz technology, where slightly absorbing media are of interest, such as for non-destructive evaluation of polymer construction materials. This technique will aid the study of internal sample structure.