1. INTRODUCTION

Terahertz (${10}^{12}\mathrm{Hz}$) light has been intensively studied as a new and promising imaging and spectroscopic analysis tool over the recent few decades. Potential applications of terahertz-pulsed time-domain spectroscopy (THz-TDS) investigated include detection of industrial defects [1], security screening [2], chemical identification [3], and cancer diagnosis [4–7]. Most THz-TDS systems are built in transmission geometry, where the detector collects the terahertz signal transmitted through the sample. However, transmission geometry can only be used to characterize partially terahertz-transparent materials. For those samples which are opaque in the terahertz range, such as hydrated biomedical materials, reflection geometry is more favorable. Furthermore, it is arguably more practical to perform terahertz sensing and imaging by detecting the terahertz beam reflected from the sample surface [8] rather than preparing a thin slice of sample to be measured in transmission.

THz-TDS in transmission geometry is a well-developed technique. Calculation of the optical properties in this geometry is highly mature and reliable. In comparison, there are still some challenges in the metrology for terahertz reflection spectroscopy. One of the main barriers for accurately extracting optical parameters in reflection geometry is the great sensitivity to any phase misalignment between the sample and reference measurements. Previously, we have demonstrated the net equations for calculating the absorption coefficient and refractive index for the reflection geometry {Eqs. (8) and (9) in [9]}. By assuming a constant reflection coefficient, we were able to calculate the change in the absorption and refractive index when a 1° phase misalignment occurs for a given phase difference. Although the absorption coefficient and refractive index are affected by both the real and imaginary part of the reflection ratio, from Fig. 2 in [9] we can see that the absorption coefficient is much more sensitive to the phase-misalignment problem: even a 1° phase misalignment will cause a significant error in the absorption coefficient calculation [10]. The absorption coefficient is an important parameter in various terahertz biomedical applications. For instance, the nonlinear trend of the concentration-dependent absorption coefficient is the key in studying the collective vibrational modes and probing the hydration shell around proteins [11,12]. As it is closely related to the hydration level of the tissue, the absorption coefficient also provides the contrast between different biological tissues such as tumor and healthy tissue [6] and fat and fibrous tissue [7]. Therefore, overcoming the challenge of the phase misalignment and accurately extracting the optical parameters in terahertz reflection measurements is of great significance in the development of terahertz applications.

The phase misalignment can be due to several reasons, such as any slight difference in position between the sample and the reference measurements, system fluctuations caused by fiber drift (if fibers are used), and mechanical jittering of the optical delay stage. Therefore, it is most severe in fiber-based terahertz systems. Several approaches have been devised to try to avoid the phase-misalignment problem, including ellipsometry [13]. However, a very high-quality polarizer is required for ellipsometry and it is difficult to obtain a polarizer with sufficient extinction ratio and transmission at terahertz frequencies. Previous works by Pashkin et al. [10], Jeon and Grischkowsky [14], Nashima et al. [15], and Hashimshony et al. [16] try to address the phase-misalignment problem by only considering the shift in sample position. Furthermore, they all have special requirements either on the sample (the sample surface needs to be optically flat) or setup (extra optical components are required). Vartiainen and co-workers employed a numerical method to solve the phase-misalignment problem by using maximum entropy [17,18]. This approach simplified the experimental procedure dramatically but was mathematically difficult.

In a typical terahertz reflection system, the sample is usually placed on top of an imaging window (Fig. 1). In many algorithms, the reflected signal from the window–air surface is taken as the reference. However, for terahertz imaging over an area of sample, the measurement time is longer and the reference is usually taken using an area near to the sample or even at a single point such as at the center point of the window. This saves measurement time but introduces noise from system fluctuations. In 1995, Thrane et al. [19] recorded the reference and sample signal in a single measurement by measuring water (the sample) in contact with a silicon window. The reflection from the air–silicon interface was used as the reference such that any phase misalignment between the sample and reference was eliminated. Jepsen et al. [20] explicitly analyzed this self-reference method in 2007 and applied this technique to determine the concentration of alcohol and sugar in commercial liquors and beverages. In their work, they first generalized this approach to be applicable to all incident angles and polarization states. A “calibration factor” was introduced to compensate for the slight pulse shape difference induced by the different beam paths between the reference and sample when the incident angle is not normal. This method is time-saving and also removes the phase-misalignment problem. However, their method requires an accurate knowledge of the imaging window thickness and has limitations when applied to imaging over sample areas as the window thickness at each spatial point can be different, as pointed out in [21]. In fact, we found that the thickness of the imaging window used in our system varied significantly with position and could not be treated as homogenous. Figure 2 shows the image of the quartz window plotted using (a) the minimum value of the signal and (b) the thickness of the imaging window calculated using the time delay between the reflections from the top and bottom surfaces and the refractive index of the imaging window. From the figures, we can see that the signal reflected from the window is nonuniform both in the amplitude and phase. The window material in our system is z-cut quartz, which is almost lossless in the terahertz range. Therefore, small differences in the thickness of the quartz window will not significantly affect the amplitude of the reflected signal in theory, but they will affect the arrival time of $E_2$.

 

Fig. 1. (a) Imaging window and (b) reflected terahertz signal.

Download Full Size | PDF

 

Fig. 2. Terahertz image of the quartz window. (a) Image of the minimum value of the reflected signal; (b) image of the window thickness calculated by the time delay between the reflections from bottom and upper surface and the refractive index of the window material.

Download Full Size | PDF

The variations illustrated in Fig. 2(a) are therefore due to the system fluctuations caused by the fiber drift, whereas the pattern in Fig. 2(b) mainly reflects the thickness variation of the quartz window itself (approximately 30 μm for a total thickness of around 2.7 mm). The thickness difference at each point of the window material will cause phase misalignment between the sample and reference data if the reference data are not taken at exactly the same points as the sample. Thus, choosing an area near the sample or any other single point as the reference is likely to cause major errors in calculating the sample optical properties. However, measuring the reference at all the same points as the sample (we refer to this algorithm as the “normal” algorithm) is obviously time-consuming and the signal fluctuation with time is not accounted for. Obviously the phase-misalignment problem caused by the window thickness could be solved by using a highly parallel window with very minimal thickness variation. However, window material of such high quality is very expensive. Therefore, in this work we provide a solution to this problem, using software to compensate for both random errors due to the system fluctuations and the systematic errors due to window thickness variation so as to accurately calculate the optical properties of the sample from terahertz reflection imaging. Our method makes use of the reflection from the air–window interface (i.e., the bottom surface of the window) and the properties of the window itself to calculate the reference. We build on Jepsen’s self-reference method to also compensate for amplitude variation due to window thickness inhomogeneity without the need to know the thickness of the imaging window at every spatial point. Our method will also simplify the experimental procedure because once the initial characterization of the window is done, no reference needs to be measured in subsequent experiments; we can simply gather all the information we need from the reflection off the bottom surface of the window during the sample measurement. In the next sections, we will introduce the theory of our proposed algorithm and demonstrate the improved accuracy achieved in extracting the refractive index and absorption coefficient.

2. THEORY

A. Baseline

We have previously reported that the reflection from the lower surface of the window material has an enduring signal (or “ringing”) which interferes with the reflection from the upper surface (the sample plane) [22]. Simply ignoring this reflection produces errors in the calculation results. The “ringing” together with the reflection from the lower surface of the window does not carry any sample information and is what we call “unwanted reflections”. After baseline subtraction, the unwanted reflections together with the fluctuations caused by system and antenna defects will be removed. The baseline is calculated from the standard complex refractive index of water (measured with terahertz transmission geometry) and a reflection measurement of water and the quartz window. All the signals used to calculate the baseline are raw signals without any windowing or filtering to preserve the original waveform of the “ringing” after the first reflection. The equations for calculating the baseline are detailed in our previous work [22]. Thus, to calculate the baseline in the “normal” method, both a reflected air signal and reflected water signal are needed at the same position as the sample measurement.

B. Transfer Function and Reflection Quotient

As illustrated in Fig. 1, there are two reflections in our signal. For the reference (air), one pulse is from the bottom of the window and the other is from the top surface. What should be noted is that the relationship between these two reflections is fixed, since it is only dependent on the properties of the window material (refractive index and the window thickness at a given position). We name this relationship the “transfer function from window to air” as shown in Fig. 3. The positional dependence is denoted by x and y, and the frequency dependence is denoted by ω. Similarly, when measuring water, the relationship between the second pulse of the air signal and the second pulse of the water signal is a function of both the properties of the imaging window and the complex refractive indices of water and air; it also will not be affected by the incident terahertz wave. We define this relationship as the “reflection quotient between water and air”. As long as the imaging window remains the same, the transfer function from window to air will not change between measurements. The reflection quotient between water and air will also stay the same. Therefore, we only need to characterize the window material and save the transfer function from window to air at every point on the imaging plane and the reflection quotient between water and air in our database once. This information can then be used to recover the reference signal and the water signal at an arbitrary point in our future measurements.

 

Fig. 3. Relationships of the transfer function from window to air and the reflection quotient between water and air.

Download Full Size | PDF

The frequency-dependent transfer function from window to air $S_{\text{window}}^{\text{air}}$ and the reflection quotient between water and air $Q_{\text{air}}^{\text{water}}$ can be calculated with

 

Fig. 4. Calibration of the THz system imaging window and calculation of reflection quotient. This step only needs to be done once.

Download Full Size | PDF

A single reference position needs to be defined to map the measurement of the transfer and quotient functions. In most practical situations, this will be provided by the origin of the mechanical stages, but in general a single fiducial mark could also do this. The transfer function from window to air is system independent and will not change over time or due to fiber drift or mechanical jittering. Therefore, this is stored in the database and can be used in all future experiments with this system. By recording the transfer function and reflection quotient, we are able to recover the air and water signal at every point of the imaging plane by simply using the first pulse of the sample signal. The first pulse should be exactly the same for both the sample and the reference for an accurate measurement but, in practice, the fluctuations in the terahertz system with time mean this is not the case in practice. Our calculated reference and water signals therefore improve the accuracy as they simulate the ideal case in which there are no fluctuations in the incident pulse. Because the reference and water signals are recovered individually for each spatial point, we avoid the phase-misalignment problem due to both the window thickness variation and the system fluctuations. The recovered reference and water signal at each sample point are calculated using

 

Fig. 5. Calculation of the corrected sample and reference signals.

Download Full Size | PDF

3. ALGORITHM VALIDATION

We demonstrate the effectiveness of our algorithm by applying it to imaging a sample of 20% isopropanol–water solution. To show that our algorithm is valid regardless of the fluctuations in the amplitude and phase of the incident wave, we repeated the sample imaging for three days. Water and the imaging window were imaged on the first day to calculate the transfer function and the reflection quotient. To prove that our algorithm can still accurately extract the optical parameters even when the amplitude of the incident wave is significantly changed, we also imaged the sample with a smaller incident wave (by lowering the voltage bias across the THz emitter). The absorption coefficient and refractive index were calculated for the samples in each day with the transfer function and reflection quotient calculated in the first day. The imaging area was $20\mathrm{mm}\times 20\mathrm{mm}$ with 0.5 mm resolution both in the $x$ and $y$ axis. Each image took 8 minutes to acquire. Each measurement was repeated 10 times in order to get a higher accuracy.

A commercial THz-TDS system from MENLO systems was used in 30° reflection geometry. The window material was z-cut quartz.

4. RESULTS AND DISCUSSION

The algorithm is implemented by a recursive function so that every spatial point on the sample goes through the processing procedure. In this way, the window thickness variation problem is taken into consideration. In this section, we evaluate our algorithm from three aspects: spatial variation, accuracy, and stability. Figure 6 is a plot of the refractive index and absorption coefficient over the imaging area at 0.45 THz. The sample is a homogeneous mixture of water and isopropanol. With the normal method, even though the sample and the reference are measured at the same spatial point, we still get obvious variations in our results as shown in Figs. 6(a) and 6(b). These variations are likely to be caused by fluctuations of the incident pulse. Our proposed algorithm, as can be seen from Figs. 6(c) and 6(d), significantly reduces this problem. The standard deviation for the refractive index is reduced from 0.03 to 0.006, giving a percentage variation of around $0.006/1.9\times 100\%\sim 0.3\%$. Similarly, the percentage error for the absorption coefficient is reduced to $\sim 0.7\%$ ($0.55/80\times 100\%$).

 

Fig. 6. Optical properties calculated with (a) and (b) the normal algorithm and (c) and (d) the proposed algorithm. The standard deviation ($\sigma$) of the sample result for each figure is written on the bottom left. The same color scales are used for (a) and (c), and (b) and (d).

Download Full Size | PDF

To test the stability of our proposed algorithm, we measured the same mixture on three different days. Since the signal variation over these three days was not large, to demonstrate that the transfer function and the reflection quotient can give us stable results with a large signal fluctuation, we also measured the mixture with a smaller incident THz pulse (obtained by reducing the bias voltage across the emitter). To illustrate the differences in the THz signal for these cases, the reflected signals from the bottom surface of the quartz window are plotted in Fig. 7. We can see that there are both phase and amplitude fluctuations. Using the transfer function and reflection quotient calculated from day 1, we extract the absorption coefficient and refractive index for the 20% isopropanol–water mixture, as shown in Fig. 8. The error bars represent the standard deviation of 676 points on the sample. The transmission measurement of the same mixture is also plotted in Fig. 8 to demonstrate the accuracy of our proposed algorithm. Even under the condition where the incident power significantly drops, as illustrated by the blue curve in Fig. 7, we can still use the same transfer function and reflection quotient to recover the reference and water signal and get a stable result of the optical properties.

 

Fig. 7. Reflections from the bottom surface of the quartz window on different days and after reducing the bias voltage applied to the THz emitter.

Download Full Size | PDF

 

Fig. 8. (a) Absorption coefficient and (b) refractive index measured by the reflection system on day 1 and calculated with the proposed algorithm in comparison with the results measured by the transmission system. The low bias curve demonstrates that even if the incident pulse has a dramatic change, the proposed algorithm can still calculate an accurate result without remeasuring the reference. (c) and (d) are the results on day 1 calculated with the normal method in comparison with the transmission measurement.

Download Full Size | PDF

To demonstrate the stability and accuracy of our algorithm, we compare our results from the reflection measurements with a transmission measurement as a benchmark. The absorption coefficient and the refractive index at 1 THz are plotted in Fig. 9. The dashed line shows the value from the transmission measurement, with the shaded area being the standard deviation of the transmission measurement. The absorption coefficient calculated from the transmission measurement is $157.88\pm 7.89{\mathrm{cm}}^{-1}$ and the refractive index is calculated as $1.94\pm 0.05$. Compared to these values, the maximum difference between the reflection measurements and the transmission measurement is 2.42% and 1.36% for the absorption coefficient and refractive index, respectively. As can be seen from Figs. 8(a) and 8(b), the result from the reflection measurements are all within the standard deviation of the transmission results. For the reflection measurements, from day 1 to day 3, the results are consistent. The low bias result is also in good correspondence with the transmission measurement. These results demonstrate that our algorithm is very accurate and stable, and even able to accommodate large system fluctuations. In contrast, as illustrated in Figs. 8(c) and 8(d), the results calculated from the normal method not only deviate from the transmission measurement, but also suffer from a larger variation on different days. It should be noted that during the validation experiments, the window was fixed with four screws at the corners to ensure that the transfer function and reflection quotient were the same in the three days. If the window must be removed (for example, for cleaning purposes), it is very likely that the transfer function and reflection quotient measured before are no longer valid. Mispositioning of the window, such as tilting, will change the optical path in the window. This will induce a $\frac{\omega }{c}\mathrm{\Delta }l$ phase misalignment to the transfer function and reflection quotient, where $\mathrm{\Delta }l$ is the change in the optical path, $c$ is the speed of light in vacuum, and $\omega$ is the angular velocity. At 1 THz, the phase misalignment caused by around 0.83 μm optical path change can be as large as 1°. Therefore, the database of the transfer functions and reflection quotients should be updated if the window is reinstalled.

 

Fig. 9. Absorption coefficient and refractive index at 1 THz calculated with the reflection measurements. The black squares are the results calculated with our proposed method, while the red triangles are calculated with the normal method. The error bars are the standard deviations. The dashed line is the value calculated with the transmission measurement, with the shaded region being the standard deviation.

Download Full Size | PDF

5. CONCLUSION

To summarize, our proposed algorithm can be used in terahertz reflection imaging systems to extract sample optical properties more accurately and efficiently than existing methods. It compensates for the phase-misalignment problem caused by the thickness variation of the imaging window, pulse drift in fiber-based systems, as well as mechanical jitters and even emitter power loss. By building on the self-referencing method to introduce a transfer function, we do not need to calculate the imaging window thickness or compensate for the beam path difference between the reference and sample. Moreover, our algorithm is applicable to any typical terahertz reflection system without any need for further modification to the system. We demonstrated the effectiveness of our algorithm by imaging a mixture of 20% isopropanol–water solution. The absorption coefficient and refractive index are calculated for each point on the imaging plane. Compared with the “normal” algorithm, our proposed algorithm reduced the standard deviation of the refractive index by a factor of 5 (from $\sigma =0.03$ to $\sigma =0.006$) and reduced the standard deviation of the absorption coefficient by a factor of 3 (from $\sigma =1.77$ to $\sigma =0.55$), resulting in percentage variation across the sample area of less than 1%. Furthermore, the result from our algorithm is highly repeatable and more accurate than the normal approach, even under the circumstance that the incident pulse significantly fluctuates in amplitude and shifts along the time axis. The maximum difference between the results from our proposed algorithm and from the transmission measurement, which is considered to be more reliable, is less than 3% at 1 THz. This accuracy is much better than existing methods and will be very beneficial to those employing terahertz reflection imaging for any application.