Invariant embedding technique for medium permittivity profile reconstruction using terahertz time-domain spectroscopy

Abstract. Reconstructing the dielectric permittivity profile (depth dependence of sample dielectric permittivity) is an important inverse problem. We present a new method for permittivity profile reconstruction based on terahertz time-domain spectroscopy signal processing. 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. Samples with known permittivity profiles have been studied experimentally using this procedure in order to verify this algorithm. This algorithm is stable to additive Gaussian white noise as shown using mathematical modeling based on the finite-difference time-domain technique. Possible applications of this permittivity profile reconstruction technique are discussed.


Introduction
Terahertz (THz) radiation is located between the infrared and microwave regions of the electromagnetic spectrum, from about 0.1 to 10.0 THz. 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. THz radiation is nonionizing in nature. Because of these properties, there are many potential applications of THz technology. 1 For instance, THz imaging and spectroscopy systems are used in the detection of concealed weapons, drugs, and explosives, 2,3 medical diagnostics, [4][5][6] and nondestructive evaluation of construction materials. 7,8 THz time-domain spectroscopy (TDS) utilizes short pulses of THz radiation with a wideband spectrum (from 0.1 to 3.0 THz) to measure the THz optical properties of various media. The electric field of a THz pulse is detected with a very high time resolution (about 50.0 fs) after the transmission of the pulse through the sample or the reflection of the pulse from the surface. By applying a fast Fourier transformation to the detected time-domain signals, it is possible to analyze the complex transmission or reflection coefficients of the sample in a wide frequency range. Time-domain data 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. 9 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.
Many ways have been presented to solve the inverse scattering problem. 9,10 However, several aspects of the inverse problem have not previously been considered, such as sample pulse function low-frequency interpolation, and sample pulse function filtering during the deconvolution process. The reconstruction results significantly depend on the particular method of solution for these problems. In addition, the stability of different algorithms with respect to noise has not been previously studied. The present work does not critically analyze previous solutions to the inverse problem. Instead, in this paper, an alternative way to solve the inverse scattering problem is presented. The results of implementing this new algorithm and its stability with respect to noise are shown in detail.
Several assumptions concerning media properties are made in the present work. The medium of interest is assumed to be nondispersive and nonabsorbing. 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 nonlinearity. When the THz penetration depth is much larger than the full optical width, a medium is considered to be low-absorbing. The dispersion of optical properties can be neglected when the variation of the medium refractive index has small fluctuations with frequency. 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). All media considered in this paper are assumed to satisfy these assumptions.
Restricting our analysis to nondispersive, nonabsorbing media, we apply the present algorithm for nondestructive evaluation of polymer construction materials. Future modifications of this method would help generalize this technique to absorbing media and media with dispersion. For example, dispersion can be described with a spectral dielectric permittivity model, such as the Debye, Lawrence, or Drude models. This improvement would aid in the development of new methods of medical diagnosis with THz TDS systems (diagnosis of skin burns, 11 study of skin structure, 12 study of tooth demineralization, 13 and cornea tissue 14 ) and methods of art inspection. 15 Most inverse scattering algorithms used in different electromagnetic spectral ranges, such as in the radio frequency range, consist of two steps for sample permittivity profile reconstruction. First, a sample pulse function RðtÞ is reconstructed based on time-domain signals. Second, the sample permittivity profile εðzÞ is reconstructed utilizing the pulse function RðtÞ. In this paper, an invariant embedding technique is suggested as the basis for the permittivity profile reconstruction procedure. This technique was previously used outside of the THz range. 16,17 Several problems should be solved in order to generalize this technique to THz TDS. These problems include the lowand high-frequency noise filtering of the sample pulse function, sample pulse function interpolation at low frequencies, and the correction of the dielectric permittivity profile using a priori information. Permittivity profile reconstruction based on an invariant embedding technique was implemented and tested experimentally and using mathematical modeling. 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.

THz Time-Domain Spectrometer
First, a THz TDS system used to collect experimental data is described. Figure 1 presents a THz time-domain spectrometer working in reflection mode. Generation and detection of THz pulses is made using optical femtosecond pulses of a fiber laser with full width at half maximum of about 80.0 fs. The generation of THz pulses is produced in a photoconductive (PC) antenna. The laser repetition rate is 50.0 MHz and the THz pulses are modulated at 102.4 kHz using PC voltage modulation. 18 The THz beam is focused on the sample. Moreover, the aperture ratio of lens is rather small: D∕f 0 ¼ 1∶2.5. Rays in THz beam transmitted through the lens have different incident angles. The maximum incident angle is equal to the aperture angle of the lens. Since the aperture angle is rather small, then the incident angle is also small. This condition is an important one for correct implementation of the algorithm. THz radiation reflects from the surface of the sample during sample signal measurements or from the reference surface during reference signal registration and is incident on the electrooptical THz field detector based on ZnTe crystal. Using a probe femtosecond optical pulse transmitted through the delay stage, the THz electrical field magnitude is detected with very high time resolution. A lock-in amplifier is synchronized with the PC antenna modulating voltage to provide a high signal-to-noise ratio.
Two signals need to be obtained for sample permittivity profile reconstruction: E s ðtÞ is the signal reflected from the sample of interest, and E r ðtÞ 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. As the method of signal detection is known, all stages of dielectric permittivity profile reconstruction can be considered.

Dielectric Permittivity Profile Reconstruction Algorithm
As previously noted, sample dielectric permittivity profile reconstruction is accomplished in two steps: reconstruction of the sample pulse function RðtÞ based on TDS system signals E s ðtÞ and E r ðtÞ and the reconstruction of the sample permittivity profile εðzÞ based on a sample pulse function RðtÞ.

Sample pulse function reconstruction
The sample pulse function RðtÞ is the medium response to the influence of an incident electromagnetic wave with a delta function amplitude time dependence. The sample Fig. 1 Schematic representation of the TDS system. pulse function is also the result of an inverse Fourier transformation applied to the sample amplitude reflectivityRðν t Þ.
Pulse function reconstruction is a signal deconvolution process. The first problem in signal deconvolution is that simply dividing the sample signal Fourier spectrumẼ s ðν t Þ by the base signal Fourier spectrumẼ r ðν t Þ to reconstruct the sample complex amplitude reflectivityRðν t Þ does not work. Since both the sample and reference signals contain no useful information at low frequencies (<0.05-0.1 THz) and high frequencies (>2.0-3.0 THz), Eq. (1) only results in noise in these spectral regions. Any sample reflectivity reconstruction procedure must contain some low-and high-frequency noise-filtering operation. The filtering procedure should allow us to obtain as much data on the sample reflectivity as possible, and should contain smoothed filter edges to prohibit the emergence of Gibbs noise in the sample pulse function. With a filter, Eq. (1) is modified as A technique utilizing a Wiener filterFðν t Þ (see Ref. 19) was applied in previous research 20,21 to process TDS system data. The formulation of a derived Wiener filterFðν t Þ with several modifications can be described as follows: (3) where Nðν t Þ is a noise power spectrum model and Sðν t Þ is a signal power spectrum model. Equation (3) is negligible in a frequency range where the noise power spectrum model Nðν t Þ exceeds the signal power spectrum model Sðν t Þ.
Conversely, if Sðν t Þ exceeds Nðν t Þ, the functionFðν t Þ is close to unity. Smooth filter edges are also provided. For simplicity, we assume a Gaussian white noise model for the noise power spectrum Nðν t Þ: whereK is a constant such that 0.0 <K < 1.0. The signal power spectrum Sðν t Þ is based on a Gaussian monopulse function f GMp ðtÞ: where ν tC is the most powerful harmonic in the monopulse spectrumf GMp ðν t Þ. The value ν tC depends on the methods that are used for generation and detection of THz pulses. The signal power spectrum model can be obtained using the following equation: where Ϝ t is the direct Fourier transform operator. Choosing a central monopulse frequency ν tC and noise spectral powerK in Eqs. (3)-(6), it is possible to provide a qualitative reconstruction of the sample reflectivity for different conditions of waveform E r ðtÞ and E s ðtÞ registration, such as the type of THz emitter and detector, the number of signal averages, and other factors. Figure 2 demonstrates an example of a Wiener filter functionFðν t Þ generated for a typical reference signal E r ðtÞ with filter parameters: ν tC ¼ 0.8 THz,K ¼ 0.002. Luckily for our TDS system, the center frequency ν tC matches a window in the atmospheric water vapor absorption spectrum. Outside this window the filtering procedure obviously leads to noise suppression in frequency ranges corresponding to high atmospheric absorption.
The second problem in the deconvolution procedure is pulse function 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 reflectivityRðν t Þ based on known frequency domain information within the ranges ½−ν t max ; −ν t min and ½ν t min ; ν t max , and to obtain lowfrequency data from the interpolation results. Since the pulse function spectrumRðν t Þ character is determined both by the internal structure of the sample and by the dispersion of media optical properties, interpolation can be easily accomplished if the investigated medium has a negligible dispersion of optical characteristics.
The modulus jRðν t Þj and phase φfRðν t Þg of the sample complex reflectivity can be interpolated separately. It is convenient to use a trigonometric interpolation series to reconstruct the amplitude reflectivity modulus jRðν t Þj. Since the sample pulse function RðtÞ is a real function, the series has the form: where Δν t ¼ ν t max − ν t min , the size of the determined sample reflectivity spectrum, with boundaries ν t min ¼ 0.1 THz and ν t max ¼ 2.0 THz. a jRj 0 and a jRj n are the interpolation coefficients: A trigonometric series is also used to interpolate the phase of the complex amplitude reflectivity φfRðν t Þg. φfRðν t Þg is an odd function, and hence, the interpolating series has the following form: The interpolation coefficients a φfRg 0 and a φfRg n can be written as The first order moment a φfRg 0 should be used in Eq. (9) instead of the zero order moment, because the sample reflectivity phase φfRðν t Þg always has a linearly decreasing character.
After reflectivity interpolation, the sample pulse function RðtÞ can be obtained using the following expressions: The interpolation type applied in this work provides good results if one deals with sharp sample permittivity profile variations or with a layered structure. Since the pulse function is known, RðtÞ, sample permittivity profile, εðzÞ, reconstruction is possible.

Sample permittivity profile reconstruction
Now consider a plane wave incident normally on the sample surface (Fig. 3). A wave reflected from εðzÞ can be described by the following integral equation: 16 where E i ðz 0 ; t 0 Þ is the incident plane wave amplitude registered at the point z ¼ z 0 , E r ðz 0 ; tÞ is the reflected plane wave amplitude registered at the point z ¼ z 0 , R þ ðz 0 ; tÞ is the scattering kernel of the integral transformation for signals registered at the point z ¼ z 0 , and z is the depth coordinate. This linear representation of the reflection process is valid when the reflected signal has no components, arising from the multiple reflections in the object layers. The kernel of the integral transformation R þ ðz 0 ; tÞ describes the region of the investigated medium in the depth range ½z 0 ; L, and depends only on the dielectric permittivity εðzÞ and effective conductivity σðzÞ of the sample. The first surface of the object is assumed to be the origin of the spatial coordinate OZ. Signals E i ðz 0 ; t 0 Þ, E r ðz 0 ; tÞ, and the integral kernel R þ ðz 0 ; tÞ can be determined if z 0 < 0, because the field amplitude can be detected only outside the object of interest. The dielectric permittivity is considered to remain constant ε 1 before the medium z < 0 and take a different constant value ε L from a certain medium depth z > L (Fig. 3). The following kernel of the integral transformation corresponds to the sample pulse function RðtÞ, obtained at the signal deconvolution step [Eq. (2)]: R þ ðz 0 ¼ 0; tÞ ¼ RðtÞ: This kernel characterizes the entire sample. For convenience, it is better to use nondimensional spatial and temporal coordinates and to normalize the integral kernels: 16 Rðx; sÞ ¼ lR þ ðz; tÞ; where l is the time needed for the wavefront to travel through the investigated medium (from 0 to L), x is the normalized optical depth (0 < x < 1), s is the normalized temporal coordinate (0 < s < 2), and cðzÞ is the dependence of the speed of light on depth. All normalized medium kernels Rðx; sÞ should satisfy the following nonlinear integro-differential equation: 16 where AðxÞ and BðxÞ are coefficients that depend on medium properties εðzÞ and σðzÞ: Equations (15) and (16) follow from dissipative wave equation.
To determine sample dielectric permittivity εðzÞ, the specified initial boundary problem in Eqs. (14), (16), and (17) needs to be solved. The integro-differential Eq. (17) is first solved numerically with the initial condition in Eq. (13) to find the AðxÞ function. For numerical computation, it is convenient to use the following notation of Eq. (17): AðxÞ: As AðxÞ was determined, the dielectric permittivity profile can be found using an inverse transformation of Eq. (16): where cð0Þ is the speed of light at the point z ¼ 0. The values cð0Þ and ε 1 should be known a priori and are often simply the values for air. Prior to permittivity profile calculation, the function AðxÞ can be corrected using some a priori information regarding the sample permittivity. For example, the sample dielectric permittivity has to be constant from a certain depth z ¼ L. To fulfill this condition, a constant should be added to AðxÞ so that the integral fðxÞ ¼ ∫ x 0 Aðx 0 Þdx 0 becomes equal to a constant at x ≥ L This correction helps to significantly reduce any error in permittivity profile reconstruction. The permittivity profile reconstruction algorithm is resumed in Fig. 4. The algorithm can be modified to take media absorbance into account. At first sight, it is impossible to reconstruct both an arbitrary dielectric permittivity εðzÞ and arbitrary conductivity σðzÞ [or absorption αðzÞ] profile simultaneously, since there is only one algorithm input function RðtÞ. However, one can use some a priori information about the sample conductivity profile σðzÞ or some approximation of the conductivity distribution during reconstruction. In addition, it is possible to connect the sample conductivity σðzÞ with sample permittivity εðzÞ utilizing a mathematical model of the complex dielectric permittivityε. In this case, reconstruction of someε parameter profile will take place instead of a direct permittivity εðzÞ and conductivity σðzÞ reconstruction. Note that the function AðxÞ in Eq. (16) depends on εðzÞ and the function BðxÞ depends on both εðzÞ and σðzÞ. We also need to include a complex dielectric permittivity modelε in the description of these functions to deal with absorbing media. After parameter profile reconstruction, one can obtain both profiles of interest fromεðzÞ.

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. A list of THz materials used in this work and their optical properties are given in Table 1. The optical thicknesses of the test samples were about 5 mm. The THz wave penetration depth in each of the materials listed in Table 1 is much higher than the sample thickness. As a consequence, material absorption can be neglected. The materials selected for the study have small dispersion for the most part of spectral range. Some of them have dispersion at frequencies above 2.0 THz. To neglect it we should reconstruct sample reflectivityRðν t Þ only in frequency range under 2.0 THz. 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. The steps for algorithm implementation for the quartz flat in air are as follows. An example of sample E s ðtÞ and reference E r ðtÞ waveforms is given in Fig. 5. The reference signal was reflected by a planar gold mirror and the sample amplitude was reflected by a 1.0-mm thick crystalline quartz window. Therefore, two pulses reflected from air/quartz and quartz/air interfaces are present in E s ðtÞ. There are no satellite pulses, which could exist due to etalon effects that arise from multiple reflections within the sample, in signal E s ðtÞ. Appearance of such pulses would lead to distortion of reconstruction results. The noise standard deviation was less than 0.5% relative to the reference peak pulse amplitude. The dynamic range and signal-to-noise ratio (SNR) 22 of reference time-domain data were 2000 and 25, respectively. Figure 6(a) illustrates the absolute value of the sample reflectivity jRðν t Þj, obtained after signal deconvolution and interpolation. Function jRðν t Þj has a strongly modulated character due to the presence of two waves in the sample waveform E s ðtÞ. The complexity of reflectivity modulation increases with an increasing number of sample pulses. Such an oscillating function can be represented as a trigonometric series during the low-frequency interpolation procedure in Eqs. (7)- (11). The sample reflectivity is distorted with noise as caused by both detector noise and changes in air composition along the THz beam path during the measurement process. Particularly strong noise lines are due to the fluctuation of water vapor content along the THz beam path.
The quartz window sample pulse function RðtÞ is presented in Fig. 6(b). The pulse function contains information about the sample interfaces as also for the sample waveform E s ðtÞ. Obviously, the pulse function RðtÞ contains two types of powerful distortions. First, there is the RðtÞ distortion caused by a sample reflectivity low-frequency interpolation error. Second, there is a Gibbs noise, which is caused by the rather sharp edges of the Wiener filter. With increasing smoothness of the filter, Gibbs noise decreases but the depth resolution of the permittivity profile reconstruction also correspondingly decreases.
The quartz plate permittivity profile reconstruction is shown in Fig. 7(a). Figure 7 also contains the results for the other test samples: Fig. 7(b) corresponds to the 1.0-mm thick quartz plate and a thick, high-resistivity floating zone silicon plate; Fig. 7(c) corresponds to the sample made of a 1.0-mm thick layer of quartz and a high-density polyethylene layer; and Fig. 7(d) corresponds to the sample made of the same quartz window and a thick polymethylpentene plate. In addition to the algorithm implementation results, any a priori information concerning sample permittivity (ideal profiles) is presented for all four figures.
The results for algorithm implementation [ Fig. 7(a)-7(d)] and the a priori data for sample permittivity allow us to conclude that profile reconstruction was produced accurately enough for all samples of interest. All reconstructed profiles contain both low-and high-frequency noise, but lowfrequency noise caused by interpolation errors leads to a larger distortion of the permittivity profile. Any small Gibbs noise present in the reconstruction results in a much smaller distortion.

Algorithm Stability in the Presence of Noise
The stability of the permittivity profile reconstruction algorithm was considered by mathematical modeling of algorithm implementation in the presence of additive white Gaussian noise in TDS system signals (Fig. 8). Models of layer medium permittivity profiles as well as models of permittivity profiles with slightly smoothed edges were considered. The interaction of THz light with matter was examined 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 εðzÞ model. First, a model for the sample permittivity profile was generated and the interaction of THz pulses with the model was considered using a finite-difference time-domain technique for solving Maxwell's equations (FDTD). 23 The assumptions made during numerical modeling were: • All media of interest were completely nondispersive and nonabsorptive. • No nonlinearities of media optical properties were present. • The structure and properties of the media did not change in the lateral directions OX and OY, and vary only with depth direction OZ. Therefore, the light interaction with εðzÞ (incident and reflected waves)  can be described with plane waves running along the OZ axis. • The time dependence of the electric field of the incident wave E inc ðtÞ was described with a Gaussian monopulse function [Eq. (5)]. • There are no satellite pulses caused by multiple reflections in the layers of the object present in the reflected signals.
After FDTD modeling, the incident E inc ðtÞ and reflected E refl ðtÞ waveforms were obtained for each studied profile εðzÞ. White Gaussian additive noise with a standard deviation σ η was added to both incident E inc ðtÞ and reflected E refl ðtÞ waveforms. Moreover, different noise standard deviations were implemented for each type of permittivity profile model εðzÞ. A permittivity profile reconstruction procedure was applied to the waveforms with noise added for estimating the permittivity profile ε 0 ðz; σ η Þ. Note, that there was no noise filtering procedure applied before permittivity reconstruction. The estimating permittivity profile ε 0 ðz; σ η Þ was compared with the initial function εðzÞ using the standard deviation of the profiles σ ε : where N is the number of permittivity profile grid points. Figure 9 contains examples of (a) the reference E inc ðtÞ, and (b) the sample E refl ðtÞ waveforms with different additive noise, as well as (c) results of permittivity profile reconstruction for various waveform noise σ η . Figure 9(c) contains both reconstruction results ε 0 ðz; σ η Þ for different noise levels and the initial permittivity profile function εðzÞ. Profiles obtained after implementation of the inverse scattering algorithm to the numerical modeling results are close to the initial profile when the additive noise standard deviation σ η is less than 1.9%. Dependence of the standard deviation σ ε for the profiles on waveform noise standard deviation σ η for the same function εðzÞ is shown in Fig. 10. 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. The quality of εðzÞ reconstruction decreases with increasing depth. Therefore, the thinner the investigated medium is, the higher the average accuracy of permittivity profile reconstruction. Some noise-filtering procedure (such as Fourier-domain or wavelet-domain noise filtering) could be implemented to the waveforms to increase the stability border to a higher standard deviation σ η (Fig. 10). Reconstruction results for all examined functions εðzÞ were stable at the same waveform noise power for noise standard deviation less than 1.9%. Certainly, the stability and uniqueness of the algorithm should be proven in analytical form, but the results shown in this chapter let us conclude that the algorithm is stable to additive Gaussian white noise.

Discussion
Experimental algorithm implementation and mathematical modeling shows high reconstruction accuracy when layered samples and samples with slightly smoothed borders between layers were studied. The average reconstruction error was smaller than 15.0% for all tested media if the waveform noise standard deviation was in the region of stability (Fig. 10). The border of stability is equal to the additive noise standard deviation σ η ¼ 1.9% and can be raised by applying a noise-filtering procedure. The theoretical limit of depth resolution is considered to be half of the TDS signal minimum wavelength, which exhibits an impulse response RðtÞ. Resolution in the range of 40 to 50 μm can be achieved for the TDS system used in this work.
The accuracy of reconstruction also depends on other factors: interpolation technique, waveform noise filtering method, media dispersion, and media absorption. The accuracy might be increased by implementing modifications to this algorithm, which may include increasing the quality of complex amplitude reflectivity interpolation, developing new methods of reconstructed permittivity profile correction, or applying a noise-filtering procedure. As mentioned above, future algorithm modifications would help to study absorbing samples and samples with dispersion of optical properties described by some complex dielectric permittivity profile model. One of the possible ways to take media absorption into account was discussed above.

Conclusions
In the present work, the ability of THz TDS applied to media dielectric permittivity profile reconstruction was shown. An algorithm for solving the inverse scattering problem based on an invariant embedding technique was developed and implemented. Problems appearing during sample pulse function reconstruction were solved. A pulse function low-and high-frequency filtering procedure based on a Wiener filter was used and trigonometric interpolation of the sample pulse function was applied. The present algorithm can be used in various fields of THz technology, where media with negligible absorption are of interest, such as for nondestructive evaluation of polymer construction materials. The results of algorithm implementation show high reconstruction accuracy for all test samples. Result of algorithm stability modeling aids in the determination of the region of TDS system noise in which reconstruction can be applied successfully. Different methods for modifying the algorithm were discussed. Future work will involve providing the ability to study absorbing media with a dispersion in optical properties. Valentin I. Alekhnovich received his PhD degree at Bauman Moscow State Technical University (BMSTU), Moscow, Russia. He is currently an assistant professor of applied mathematics chair (Fundamental Sciences Department) and an assistant professor of laser and optoelectronic systems chair (Department of Radio Electronics and Laser Technology) at BMSTU. He is the author of over 100 publications and conference presentations. His research interests include laser spectroscopy, inverse and ill-posed problems in optics, optical, and laser devices.