## Abstract

Terahertz time domain spectroscopy allows for characterization of dielectrics even in cases where the samples thickness is unknown. However, a parameter extraction over a broad frequency range with simultaneous thickness determination is time consuming using conventional algorithms due to the large number of optimization steps. In this paper we present a novel method to extract the data. By employing a three dimensional optimization algorithm the calculation effort is significantly reduced while preserving the same accuracy level as conventional approaches. The presented method is even fast enough to be used in imaging applications.

© 2011 OSA

## 1. Introduction

The progress of terahertz (THz) technology in the recent decades yields various system approaches. On the one hand side, there are continuous wave (CW) systems that operate at a single narrow band frequency [1–4]. Unfavorable here is that the periodicity of the CW signal induces ambiguous phase information which limits a reliable characterization of the samples under investigation. To overcome this it is necessary to measure the sample at multiple different frequencies, either by tuning the system and perform several measurements [5] or by utilizing a multi-color CW THz system [6].

Better suitable are THz spectrometers that exhibit a pulsed waveform [7–10]. Especially the THz time domain spectroscopy (TDS) (see [10] for a comprehensive review) is a powerful tool since the signal contains broadband information. Thus, a single measurement reveals the sample's dielectric properties over a broad frequency range. Moreover, multiple reflections from the samples facets result in so called Fabry-Pérot (FP)-echo pulses which contain additional information. Considering them during the data analysis enables for the simultaneous determination of the samples thickness and the dielectric optical parameters [11–13]. To this end, a numerical optimization procedure is employed to calculate a theoretical transfer function that describes the measurement. In detail, the procedure is the following: First a rough guessed value for the thickness is assumed and initial values for the real and imaginary part of the refractive index are deducted in the frequency range of interest. Afterwards, the theoretical transfer function is calculated individually for every single frequency step and compared to the measured transfer function. A two-dimensional numerical optimization minimizes the deviations between the calculated and the measured transfer function. This leads to the complex refractive index of the sample at a single frequency. To extract the broadband dielectric function, the optimization has to be performed successively for every frequency step.

However, the calculated material parameters are a function of the guessed thickness. If this deviates from the reality, the extracted refractive index and absorption coefficient is compromised by FP oscillations [11–13]. Only for the exact thickness the FP oscillation vanishes. Performing the optimization for a set of thicknesses and minimizing the FP oscillation allows consequentially for determining the sample's thickness and the dielectric properties simultaneously. This data extraction procedure is well suited for spectroscopic applications. Yet, the high number of optimizations (number of frequency steps times number of investigated thickness steps) makes this approach slow and not suitable for analyzing data from measurements where high speed is required like in THz imaging or in most industrial applications.

In this paper, we present a novel method to extract the data. Instead of using a large number of successive optimization procedures we employ a single more sophisticated optimization. To this end, we define the broadband complex transfer function from a three dimensional set of parameters. These are i) the sample thickness, ii) the real part, and iii) the imaginary part of the broadband refractive index. We optimize for the linear combination of these three parameters that best agrees with the measured transfer function. Since the measurement ranges over a broad frequency interval, it is necessary to define the complex refractive index also for the entire frequency window. Here, we employ initial values derived in a preprocessing step as described in the next section.

The proposed method outperforms the conventional algorithms significantly in terms of speed while the results are comparable. The required calculation time is only in the range of a few milliseconds and thus more than one thousand times faster as conventional approaches. As state of the art TDS systems have a scan rate of about 100 Hz [14], the proposed algorithm can be referred to as real-time signal processing since the evaluation time is in the same range as the signal acquisition time. Therefore, the proposed method is ideally suitable for THz imaging and other applications where high speed and precise accuracy is required.

## 2. Algorithm

A sample can be characterized in TDS transmission geometry by measuring the complex spectrum of the THz-wave propagating through the sample and a reference spectrum. The transfer function of the sample is given by the ratio between these two complex spectra [10]. From this transfer function *H*(*f*) we derive initial values for the complex refractive index in a first step. Here, we use the conventional equations that neglect FP echoes and the phase of the complex reflection coefficients [15]. For an assumed initial value of the thickness *L _{0}*, the frequency-dependent initial values for the real and imaginary part of the refractive index

*n*and

_{0}*κ*, respectively, can be approximated by:

_{0}*c*denotes the vacuum speed of light. However, the FP echoes induce a FP oscillation on these parameters. The oscillation period

_{0}*f*is given by [13]:Since this oscillation is not related to the physical material parameters, we filter it out using a band stop filter that is centered on

_{FP}*f*. This results in smooth material parameters as illustrated in Fig. 1 . Please note that these filtered values are only employed as initial values in the optimization procedure. To extract the optical parameters in the final step, the raw measured data is used.

_{FP}Yet, if the sample exhibits absorption features which occur at a frequency spacing in the range of the FP oscillations, the data extraction might be compromised. As example, alpha-lactose features absorption lines at 500, 1200, 1400 GHz, and above. But if the sample is e.g. 1 mm thick, the FP resonance has a period of only 90 GHz. Thus, these oscillations can easily be filtered out while preserving the dielectric features of the sample. However, one can imagine cases, where a thin sample exhibits absorption features close to the FP resonance. In the case of alpha-lactose a sample thickness of around 200 µm would be critical. Here, the algorithm might fail - along with all the other methods as described in [11–13]. In general, the lower limit for a robust thickness determination with THz TDS depends on the dispersion as well as the absorption features of the sample under investigation. While the lower limit for alpha-lactose can be seen around 200 µm, for weakly dispersive materials as silica or most polymers it is possible to analyze samples even in the sub 100 micron regime [13].

The derived initial values, real part *n* and imaginary part *κ* of the refractive index together with the thickness *L*, constitute the parameter set to calculate the theoretically complex transfer function that describes the sample response:

Here, *M* is the maximal FP echo pulse that appears in the time window of the measurement. For the optimization we use three scalar parameters *ξ*, *ψ* and *ζ*:

The reason for this definition is that the induced phase shift by the sample does not scale with *n* but is related to (*n*-1). Therefore, the order of magnitude of the three parameters is comparable, which is favorable for numerical optimization. The resulting theoretical transfer function is given by:

The next step is to perform a numerical optimization where the Nelder-Mead simplex algorithm [16] is used. Here, we define an error function *ΔH* as:

As the real and the imaginary part of the transfer function is considered, all three values, *n*, *κ* and *L*, significantly influence the error function. Thus, the optimization is expected to derive the best values for all three fitting coefficients *ξ*, *ψ*, and *ζ*.

After obtaining the optimal values for *ξ _{opt}*,

*ψ*, and

_{opt}*ζ*, the thickness

_{opt}*L*is known and the resulting material parameters can be calculated:

_{opt}However, the resulting optical constants *n _{opt}* and

*κ*originate from the filtered initial values. Thus, absorption features of the material under investigation may be corrupted due to the filtering. To obtain precise material parameters, we apply in the final step a single successive optimization similar to [13] using the derived thickness

_{opt}*L*. Here,

_{opt}*n*(

*f*) and

_{j}*k*(

*f*) are optimized at every frequency step

_{j}*f*using the error function:

_{j}## 3. Results

To validate the capability of the proposed method, we apply the algorithm to measured data and compare the results with conventional procedures. To this end, we measure three different samples: a silica wafer that is weakly dispersive, a slightly dispersive sample of polyethylene terephthalate (PET), and an alpha-lactose-polypropylene (PP) compound exhibiting pronounced absorption features. To record the data we employ a standard THz TDS setup driven by a titanium-sapphire laser (see e.g., [10] for a description). The repetition rate of the laser is 80 MHz and the pulse duration is 80 fs.

For the silica sample, we assume an initial value of about 750 µm. The corresponding initial material parameters are shown in Figs. 2(a) and 2(b). As can be seen in this figure, the FP oscillations of about 100 GHz tower over dispersion effects. The filtered data are also shown in Figs. 2(a) and 2(b). These values are used as initial values for the optimization procedure. The optimization leads to a thickness of 707.7 µm after a computation time of 12 ms. Using the conventional algorithm [13], we obtain an optimum thickness between 707 and 708 µm as shown in Fig. 2(c). Here, the remaining strength of the FP oscillation as shown in the figure is used as figure of merit to determine the correct thickness [13]. For the conventional method the computation time is about one minute at the same computer. The flat valley of about 1 µm results from the limited signal to noise ratio of the measurement that also affects the accuracy of the thickness determination.

The thickness obtained by the new concept agrees very well with the conventional procedure. The extracted material parameters are shown in Figs. 2(a) and 2(b). The deviation between the extracted data based on the filtered initial values and the successive optimization is almost negligible. Thus, if speed is in favor over accuracy the last step of the algorithm may be skipped for an even faster data analysis. This way the calculation time is reduced based on the number of investigated frequency steps by around a factor of four.

To give an inside to the algorithms outcome, we calculate the two dimensional error function that results when setting one of the parameters (*ξ*, *ψ*, *ζ*) constant (Figs. 3(a)
–3(c)). The figure shows that a distinctive minimum occurs for values close to the real physical ones. Particularly, the choice of *ξ* and *ζ*, which denote the thickness and the refractive index of the sample, respectively, influence the error function noticeably. Due to the definition of *ξ* and *ζ* (cf. Eq. (5)) circularly shaped sets of local minima result. Yet, the remaining error value of the global minimum that corresponds to the real physical thickness is significantly lower than the errors for the other minima. The parameter *ψ*, as measure for the imaginary part of the refractive index, only weakly influence the error function in this case due to the relatively high transparency of the sample. The global minimum is given for a defined value of *ψ* though.

For the other two samples, we chose initial thicknesses of 320 µm (PET) and 2400 µm (alpha-lactose-PP compound), respectively. The resulting initial values and extracted parameters are shown in Fig. 4 and Fig. 5 . In all cases, the results of the new method obtained in a few milliseconds agree well with the conventional ones that are computed on a close to a minute scale.

While the silica wafer exhibits an almost vanishing dispersion and very low absorption coefficient within the investigated frequency range, the material parameters of PET are noticeably frequency-dependent. Due to the selective filtering of the FP oscillations during deriving the initial values, the dispersive shape of extracted parameters is not compromised.

The last sample measured, the alpha-lactose-PP compound, exhibits absorption features at 0.5, 1.2 and 1.4 THz. As the FP resonance is only 40 GHz due to the relatively thick sample (2.6 mm) it is possible to filter it out easily without compromising the absorptive behavior of the compound.

The new method clearly delivers comparable results to the conventional approaches in fractions of seconds. Yet, convergence issues have to be discussed. For that we vary the initial values and observe the algorithms outcome. As example, we chose the PET sample due to its moderate but not vanishing dispersion. Since the initial values for *n _{0}* and

*κ*result from the choice of

_{0}*L*, we vary this parameter. We iteratively perform the computation for initial values of the thickness between 30 µm and 600 µm. The resulting thickness of the algorithm is shown in Fig. 6 together with the remaining error function. As one can see in the figure the algorithm converges to the correct value within the range between 30 and 451 µm. For other initial values, it converges to another local minimum. Yet, the remaining value of the error function allows us to distinguish between the global minimum and other minima. Using more sophisticated optimization algorithms that are more robust against local minima distortion would therefore extend the convergence range for the initial values.

_{0}While in most cases, the thickness of the sample under investigation is approximately known, there are applications where no information of the samples thickness is given. For that the algorithm could be used iteratively for different initial thickness values. The result with the lowest error function is than defined as the algorithms outcome. This concept is further discussed in [17]. Further work will be concentrated on improving the robustness of the algorithm by employing optimized error functions and optimization methods.

## 4. Conclusion

In conclusion we present a new method to analyze data obtained by terahertz time domain spectroscopy. The method is based on a three dimensional optimization procedure where a broadband error function between the measured and a theoretical calculated transfer function is minimized. Doing so allows us to simultaneously determine the dielectric properties of the samples under investigation as well as their thickness in fractions of seconds. We compare the new procedure with conventional algorithms and show a speed enhancement in the range of three orders of magnitude.

## Acknowledgments

I acknowledge Thorsten Probst and Benedikt Scherger (Philipps-Universität Marburg) for their support in acquiring the measured data as well as Martin Koch (Philipps-Universität Marburg) for fruitful discussions about the manuscript.

## References and links

**1. **S. Matsuura, M. Tani, and K. Sakai, “Generation of coherent terahertz radiation by photomixing in dipole photoconductive antennas,” Appl. Phys. Lett. **70**(5), 559–561 (1997). [CrossRef]

**2. **S. Verghese, K. A. McIntosh, S. Calawa, W. F. Dinatale, E. K. Duerr, and K. A. Molvar, “Generation and detection of coherent terahertz waves using two photomixers,” Appl. Phys. Lett. **73**(26), 3824–3826 (1998). [CrossRef]

**3. **K. J. Siebert, H. Quast, R. Leonhardt, T. Löffler, M. Thomson, T. Bauer, H. G. Roskos, and S. Czasch, “Continuous-wave all-optoelectronic terahertz imaging,” Appl. Phys. Lett. **80**(16), 3003–3005 (2002). [CrossRef]

**4. **M. Scheller, J. M. Yarborough, J. V. Moloney, M. Fallahi, M. Koch, and S. W. Koch, “Room temperature continuous wave milliwatt terahertz source,” Opt. Express **18**(26), 27112–27117 (2010). [CrossRef]

**5. **X. Wang, L. Hou, and Y. Zhang, “Continuous-wave terahertz interferometry with multiwavelength phase unwrapping,” Appl. Opt. **49**(27), 5095–5102 (2010). [CrossRef] [PubMed]

**6. **M. Scheller, K. Baaske, and M. Koch, “Multifrequency continuous wave terahertz spectroscopy for absolute thickness determination,” Appl. Phys. Lett. **96**(15), 151112 (2010). [CrossRef]

**7. **D. H. Auston and K. P. Cheung, “Coherent time-domain far-infrared spectroscopy,” J. Opt. Soc. Am. B **2**(4), 606–612 (1985). [CrossRef]

**8. **B. Sartorius, H. Roehle, H. Künzel, J. Böttcher, M. Schlak, D. Stanze, H. Venghaus, and M. Schell, “All-fiber terahertz time-domain spectrometer operating at 1.5 μm telecom wavelengths,” Opt. Express **16**(13), 9565–9570 (2008). [CrossRef] [PubMed]

**9. **M. Walther, B. M. Fischer, A. Ortner, A. Bitzer, A. Thoman, and H. Helm, “Chemical sensing and imaging with pulsed terahertz radiation,” Anal. Bioanal. Chem. **397**(3), 1009–1017 (2010). [CrossRef] [PubMed]

**10. **P. Jepsen, D. G. Cooke, and M. Koch, “Terahertz spectroscopy and imaging—modern techniques and applications,” Laser Photonics Rev. **5**(1), 124–166 (2011). [CrossRef]

**11. **L. Duvillaret, F. Garet, and J. L. Coutaz, “Highly precise determination of optical constants and sample thickness in terahertz time-domain spectroscopy,” Appl. Opt. **38**(2), 409–415 (1999). [CrossRef]

**12. **T. D. Dorney, R. G. Baraniuk, and D. M. Mittleman, “Material parameter estimation with terahertz time-domain spectroscopy,” J. Opt. Soc. Am. A **18**(7), 1562–1571 (2001). [CrossRef]

**13. **M. Scheller, C. Jansen, and M. Koch, “Analyzing sub-100µm samples with transmission terahertz time domain spectroscopy,” Opt. Commun. **282**(7), 1304–1306 (2009). [CrossRef]

**14. **I. Duling and D. Zimdars, “Terahertz imaging: revealing hidden defects,” Nat. Photonics **3**(11), 630–632 (2009). [CrossRef]

**15. **P. U. Jepsen and B. M. Fischer, “Dynamic range in terahertz time-domain transmission and reflection spectroscopy,” Opt. Lett. **30**(1), 29–31 (2005). [CrossRef] [PubMed]

**16. **J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence properties of the nelder–mead simplex method in low dimensions,” SIAM J. Optim. **9**(1), 112–147 (1998). [CrossRef]

**17. **M. Scheller and M. Koch, ““Fast and accurate thickness determination of unknown materials using terahertz time domain spectroscopy,” J. Infrared Milli. Terahz. Waves **30**(7), 762–769 (2009). [CrossRef]