An introduction to terahertz time-domain spectroscopic ellipsometry

In the past, terahertz spectroscopy has mainly been performed based on terahertz time-domain spectroscopy systems in a transmission or a window/prism-supported reflection configuration. These conventional approaches have limitations in regard to characterizing opaque solids, conductive thin films, multiple-layer structures


I. INTRODUCTION
Polarization is a fundamental physical property of light, which can be altered upon interaction with a medium. In other words, the polarization state of light contains information about the interacted sample. Polarimetry is a technique in which the properties of a sample are retrieved from changes in the polarization state of light. This can be done for either transmission or reflection at oblique incident angles, while reflection-type polarimetry is more commonly applied due to its ability to measure opaque samples. Reflection polarimetry is a branch of ellipsometry, although in many cases, they are synonymous. In this Tutorial, we will use the term ellipsometry as we will focus on sample characterization in a reflection configuration.
To correlate the polarization change with the sample properties using the method of ellipsometry, one compares the p and s reflection coefficientsrp andrs and expresses the ratio asρ =rp/rs ( ∼ represents complex values). When the incident light has equal p and s components, this ratio is correlated with the complex reflected electric fields in the p and s directions,Ẽrp andẼrs, respectively, as follows: 1ρ where tan Ψ and Δ represent the magnitude ratio and phase difference, respectively, ofρ.Ẽip andẼis are the incident electric fields in p or s directions, respectively. By establishing an optical model to describe the light-matter interaction,rp andrs, and, hence,ρ, can be predicted by Fresnel coefficients as a function of the sample properties, such asρ(ñ) =Ẽrp/Ẽrs, whereñ is the complex refractive index of the sample. In contrast, characterization is done by solving the equation either analytically, numerically, or by model fitting. Therefore, ellipsometry extracts the sample properties by self-referencing the two orthogonal electric fields from the sample, without the need for an extra reference signal.
The first use of the term "ellipsometry" dates back to 1945. 2 At the present time, spectroscopic ellipsometry has become a mature characterization technology in the infrared (IR)-UV range, with many commercial ellipsometer products available. It is straightforward to consider extending this technique further to the terahertz (THz) wavelengths, which is desirable in this regime compared to conventional THz spectroscopy methods, and it was first demonstrated by Nagashima and Hangyo in 2001 based on THz-TDS (time-domain spectroscopy). 3 Conventional THz spectroscopy uses a transmission or window/prism-based reflection configuration, which allows the study of transparent bulk materials, most liquids, and soft tissues. However, absorptive solids and conductive thin films, such as most doped semiconductors, amorphous inorganic materials, and conductors, are challenging to be precisely characterized. This is because in transmission, absorptive solids attenuate the light very quickly while conductive thin films are usually too thin to cause sufficient phase change. In reflection, these solid samples cannot contact well with a supporting medium that an uncontrollable air gap will be induced. The noncontact reflection scheme, which compares the reflection from the sample with the reflection from a reference medium, is rarely used for characterization because of the "phase uncertainty" problem 4 caused by the height difference between the sample and the reference medium. For these samples, ellipsometry provides an ideal noncontact and self-reference modality. Actually, ellipsometry is a versatile technology that is not limited only to these sample types; it is also powerful in characterizing anisotropic samples 5 and multiple-layer structures 6 and in investigating magneto-optical effect [7][8][9] and polarization-sensitive devices. [10][11][12] THz spectroscopic ellipsometry has been done in either frequency-domain (FD) or time-domain (TD). Due to large discrepancies in the instrumentation between these two techniques, they have very different characteristics in regard to the bandwidth, detected quantities, measurement methods, and data processing methods. Actually, within the class of FD ellipsometers, there are numerous source-detector combinations, such as BWO (backward wave oscillator)-Golay cell, 9,13,14 BWO-SNA (scalar network analyzer), 15 and black-body source in FTIR (Fourier transform infrared) plus synchrotron-bolometer. 16 A majority of these THz FD ellipsometers were proposed by the group in University of Nebraska-Lincoln and their co-workers. On the contrary, the category of THz-TDS ellipsometers shares similar operation mechanisms. The studies reported so far have mostly used nonlinear crystals or photoconductive antennas (PCAs) for THz emission and detection, [17][18][19][20] and very recently, the use of air-plasma filament polarimetry has also been demonstrated. 21 As a tutorial aiming at delivering the key technical guidance and skills for a specific technique, we will only focus on ellipsometry based on THz-TDS due to the huge differences from the FD systems.
THz-TDS ellipsometers have some advantages and challenges compared to THz-FD ellipsometers and commercial IR-UV ellipsometers. The first obvious merit is the coherent detection of electric fields, which simultaneously provides both the amplitude and the absolute phase in a ultrabroad bandwidth at a fast acquisition rate (e.g., >20 Hz for systems using a femtosecond laser at a repetition rate of 100 MHz). Coherent detection also simplifies the polarization control, reducing the number of measurements needed. In detail, tan Ψ and Δ can be directly obtained from THz-TDS using a polarizer-sample-analyzer (PSA) scheme with p and s reflections measured using two analyzer orientations. In contrast, most FD ellipsometers are intensity-based, using measurements of at least four polarization directions, where typically the p, s, and ±45 ○ components are needed to extract tan Ψ and cosΔ, and further measurements of the left and right circular polarization components manipulated by additional phase compensators are required to determine Δ in the range of [−180 ○ , 180 ○ ] from cos Δ. 1 Broad bandwidth is another advantage this method offers compared to THz FD sources. Typical TDS systems cover 0.1-4 THz, while most THz CW sources have a limited tunable range and require additional frequency-multipliers to achieve a broader bandwidth. Nevertheless, their excellent spectral resolution down to 1 MHz could be useful in special applications to resolve fine resonant features. 14,22 Finally, the picosecond-temporal resolution (hence, the absence of standing-wave issue 14 ) enables selecting specific reflective pulses in the time domain, simplifying the data processing, and may provide additional spectral information. A general downside to using THz ellipsometry is the beam divergence issue, which is physically unavoidable when the wavelength of the light increases to the quasi-optical regime. We will discuss this in detail in Sec. IV B. A particular weakness of THz-TDS ellipsometers is the high sensitivity of the single-pixel detector to the THz beam alignment. Most TDS detectors require the THz beam to be precisely focused; hence, multiple reflections from bulk samples could be out-of-focus and have different detection sensitivities, making the spectrum containing all these reflections differ from the theoretical reflection model. This can be solved in many cases by temporally removing the high-order reflections. Pulse shift is another special issue that could be associated with the time-domain detection modality, causing phase errors especially at high frequencies.
In this Tutorial, we will focus on these unique characteristics of THz-TDS ellipsometry and provide key technical guidance and skills to obtain accurate ellipsometric measurements. We will discuss the theory in Sec. II, fundamental optics in Sec. III, error analysis and calibration in Sec. IV, and sample characterization in Sec. V. Finally, experimental demonstration will be presented in Sec. VI.

A. Categories
Ellipsometry can be classified into two categories: standard ellipsometry and generalized ellipsometry. Standard ellipsometry refers to measurements with no cross-polarization induced, that is, no s reflection is produced under a p incidence or vice versa and no depolarization occurs. Isotropic materials and uniaxial anisotropic materials whose optical axis is specially aligned parallel or perpendicular to the incident plane can be measured by using this technique. The most commonly reported THz-TDS ellipsometric measurements are carried out via standard ellipsometry. For anisotropic materials with randomly orientated optical axis, conversion between the s and p polarizations occurs, which requires generalized ellipsometry in order to measure cross-polarization reflections. If depolarization occurs, such as by materials with a wavelength-comparable surface roughness, Mueller-matrix (introduced next in Sec. II C) generalized ellipsometry has to be applied to describe partially polarized or unpolarized APL Photonics TUTORIAL scitation.org/journal/app light. Generalized ellipsometry in the THz regime has been mostly demonstrated by FD systems, 5,14,16,23 while THz-TDS ellipsometers can also be built in this form to measure cross-polarization reflections.

B. Fresnel coefficients
Fresnel coefficients describe the reflection and transmission of light for an established optical model. Here, we only introduce reflection coefficients used in reflection-type ellipsometry. For bilayer structures containing two semi-infinite media 1 and 2, as shown in Fig. 1(a), the p and s reflection coefficients for light incident from medium 1 can be expressed as follows: rs =Ẽ rs Eis =ñ 1 cos θ 1 −ñ 2 cos θ 2 n 1 cos θ 1 +ñ 2 cos θ 2 , whereñ 1 (θ 1 ) andñ 2 (θ 2 ) are the complex refractive indices (incident/refracted angles) in medium 1 and medium 2, respectively. Snell's law connects θ 1 and θ 2 using the following relation: For an optical model containing three layers, from 1 to 3, with the middle layer 2 having a finite thickness d 2 , as shown in Fig. 1(b), the corresponding p(s) reflection coefficient for light incident from medium 1 is given by where β = 2πd 2 λñ 2 cos θ 2 (6) is called the film phase thickness.r 12,p(s) andr 23,p(s) in Eq. (5) are p or s reflection coefficients from medium 1 to 2 and from 2 to 3, respectively, calculated by Eq. (2) or Eq. (3) according to the polarization state. For a stratified medium containing more than three optical layers, we can apply Eq. (5) for every three layers iteratively. For example, for a four-layer structure containing media 1-4 with light incident from medium 1, as shown in Fig. 1(c), we can first calculate the reflection of the lower three layersr 234 by using Eq. (5). The total reflectionr 1234 is calculated again using Eq. (5) by replacing the lower reflection withr 234 : r 1234,p(s) =r 12,p(s) +r 234,p(s) exp(−2iβ) 1 +r 12,p(s)r234,p(s) exp(−2iβ) .
The same principle can be applied for an arbitrary number of layers. Using these reflection coefficients, we are able to calculate the polarization state of the reflected light as a function of the sample properties [i.e.,ρ(ñ)]. The sample properties are determined from the best fit between the calculated and measured results.

C. Data expression
An ellipsometer may contain multiple polarization-dependent optical elements. The measured quantity depends on the orientation of these components. It could be difficult to directly express the measured signal as a function of the sample properties when there are numerous polarization variations during the propagation. The Jones matrix and Mueller matrix are the two mathematical approaches used to express polarization-dependent light propagation. A Jones matrix J is a 2 × 2 complex-valued matrix that relates the input and output electric fields (expressed as Jones vectors) as follows: A Mueller matrix M is a 4 × 4 real-valued matrix. In this case, the input and output are expressed as Stokes vectors S whose parameters are intensity-based quantities. The relationship becomes APL Photonics TUTORIAL scitation.org/journal/app where the Stokes parameters are related to the intensities as follows: where the subscripts x, y, +45 ○ , −45 ○ , R (right-circular), and L (leftcircular) represent the polarization components. Optical elements, sample response and coordinate rotation, can be expressed by independent Jones matrices in the sequence of their positions in the propagating path as J 1 J 2 J 3 . . . or by Mueller matrices as M1M2M3. . .. The product of these matrices relates the input from the source and the output to the detector. The elements of the final product matrix contain the sample response, which are functions of the sample properties that can be estimated by Fresnel coefficients. The major difference between the Jones matrix and Mueller matrix is that only the Mueller matrix can express unpolarized and partially polarized light. For polarized light, both expressions are mathematically identical, while because the Mueller matrix deals with intensitybased Stokes parameters, it is sometimes more favorable for FD ellipsometry in which the intensities of different polarization components are measured. For THz-TDS ellipsometry, which involves directly measuring the complex electric fields without considering depolarization, the Jones matrix is much more convenient. Actually, the coherent detection modality of THz-TDS significantly simplifies the polarization control and the measurement (see Sec. III B). As a result, the relationship between the measured electric fields and the sample properties [i.e., Eq. (1)] can be easily derived, and matrix multiplications in terms of simple measurements becomes unnecessary. 18,[24][25][26] Nevertheless, the Jones matrix can be used to assist the analysis of more complicated situations, such as generalized ellipsometry or when taking imperfections of optical elements into account. 20,27

A. Beam control
Ellipsometry retrieves the sample properties from the observed change in polarization, mathematically based on tan Ψ and Δ. The increased measurement accuracy relies on the high sensitivity of tan Ψ and Δ to the sample properties. For bulk materials, this is achieved near the Brewster angle (or pseudo-Brewster angle for absorptive samples) θB. Assuming a medium withñ = n − iκ, Figs. 2(a) and 2(b) show the calculated tan Ψ vs the incident angle θi from air to the medium with n varied from 1.5 to 3 (κ fixed at 0) and with κ varied from 0 to 1 (n fixed at 2), respectively. The biggest variation of tan Ψ between curves at a fixed angle is found around θB, while the smallest variation is observed at the normal incident and 90 ○ . Similar characteristics can also be observed for Δ (not shown). The analysis reflects two characteristics of ellipsometry. First, the incident angle should be adjusted specifically for a certain sample according to its properties. Doing this requires a robust, ideally, electrical control of the incident angle. Second, the angular region with a high sample sensitivity is usually accompanied by a high angular sensitivity as well, especially for materials with a large refractive index, as can be seen in the rapid variation of tan Ψ at lager incident angles in Fig. 2(a). This means the incident angle should be precisely set and measured. To satisfy these requirements, fiber-coupled photoconductive antennas (PCAs) are highly recommended for THz generation and detection. Free-space coupled PCA, optical rectification, and electro-optic sampling are sensitive to the optical alignment of the pumping and probing beams; hence, the emitter and detector are usually fixed. The incident angle is changed by rotating off-axis mirrors, which in turn requires realigning the THz optics. 20,25,26 In comparison, fiber-coupled antennas can be freely moved without affecting the coupling between the femtosecond laser beams and the antennas. 18,19,[27][28][29] The optics on the emission and detection sides can be assembled on two independent rails, respectively, as shown in Figs. 3(a) and 3(b). In this way, the incident angle can be easily adjusted by rotating the rails around the focal point on the sample, without causing a significant misalignment, as indicated by the gray arrows in Fig. 3(a). As the polarization-dependence is typically weak at small incident angles, designing a rotational stage with a switchable range of 45 ○ -90 ○ is suitable for most sample types, which also allows additional transmission spectra to be measured for transparent samples. Figure 3(b) shows the optical arrangement from the incidentplane view. For PCAs used as an emitter and detector, a pair of lenses L1 and L4 (could be replaced with parabolic mirrors), typically with a f -number around 1.3-2, are set next to the antennas to collimate the beam from the emitter and to focus the collimated beam to the detector. Another pair of lenses, L2 and L3, is used to focus the beam onto the sample and to collimate the reflected beam from the sample. The selection rule for this pair of lenses will be discussed in Sec. IV B.

B. Polarization manipulation
Polarization state measurement is another key feature of an ellipsometer. We have discussed earlier that the coherent detection of THz-TDS greatly simplifies measurements by only detecting two orthogonal electric fields. Therefore, unlike IR-UV ellipsometers that require additional phase compensators to measure the circular polarization components, THz-TDS ellipsometry only needs polarizers in a PSA configuration to obtain the absolute polarization state. Typically, three polarizers are necessary for a precise polarization manipulation in standard ellipsometry, as shown in Fig. 3 Most THz PCAs are quasi-linearly polarized in emission and detection. 12 The degree of linearity strongly depends on the antenna design, but most of them still have considerable cross-talk. Therefore, two polarizers, P1 and P3, are placed in front of the emitter and detector, respectively, with the passing direction aligned parallel to the main-polarizing direction of the antennas. 28,29 As we will prove in Sec. IV C, it is convenient to regard this polarizer-antenna combination as a single unit, providing an ideal linear emission and detection. They are normally set at 45 ○ (to the s direction). Another polarizer, P2, acts as the analyzer to select the reflected s or p component. The polarization states at different beam positions are illustrated as insets in Fig. 3(b). The emission unit sends a perfectly 45 ○ linearly polarized light after P1, containing equal s or p components. Sample reflection causes a polarization variation, which is reconstructed by rotating P2 to measure the s and p components, respectively. The detector unit aligning at 45 ○ ensures that both s and p components after P2 can be detected with equal sensitivity. For generalized ellipsometry in which crosspolarization occurs, an additional polarizer can be placed on the emitter side after P1 to set the incident light to the sample as s or p polarized. In this way, four electric field quantities of Erpp,Ẽrsp,Ẽrps, andẼrss can be measured, and three relationships, typically by normalizing the former three electric fields toẼrss, can be obtained. Note that in any measurement, it is not recommended to rotate the antenna units as any mechanical adjustment to the emitter and detector will lead to a sensitive change in the optical alignment, hence the measured electric fields cannot be self-referenced.

A. Measurement error
Measurement errors occur in any practical experiment, ultimately limiting the accuracy that can be achieved. Error-propagation analysis is essentially a sensitivity analysis, which is especially important for ellipsometry as the sensitivity highly depends on the incident angle, and it is particularly useful for samples without a clear definition of θB (e.g., conductive materials). A poor sensitivity significantly magnifies the measurement error in regard to the characterization results (e.g., measured at near normal incidence). Error-propagation calculation is also important for result analysis. THz-TDS measurements offer a frequency-dependent signal-tonoise ratio (SNR). Analyzing the noise error enables the determination of the data credibility and the available bandwidth also. Table I provides steps to perform error-propagation analysis. If the analysis is performed prior to an experiment to analyze the sensitivity, the first two steps are required.
Step 1 estimates the sample propertiesñ sample (ω), which can be obtained from a proper dielectric model or literature. The accuracy requirement is low for sensitivity analysis. Based on the optical model corresponding to the measurement,rp(ω) andrs(ω) can be calculated usingñ sample (ω).
Step 3 adds random noise to the spectrum, as shown in the equations. If the analysis is performed after an experiment,Ẽrp andẼrs are directly obtained from the measurement. N(ω) in the equation represents the noise amplitude randomly assigned within [0, Nmax], where Nmax is the maximum noise amplitude that can be estimated from the noise floor. ϕ(ω) is the noise phase randomly distributed within [0, 2π]. Note that here we only consider the white noise mainly from the detector, hence, Nmax is independent of the reflectivity and frequency. In this case, a weak reflection leads to a small signal and, hence, a low SNR.ρnoise(ω) is regarded as the noise-affected ratio used to characterizeñ sample in Step 4 (characterization methods are detailed in Sec. V). Steps 3 and 4 are repeated M times (e.g., M > 100) to enable statistical analysis.
Based on these steps, we show an example of error analysis for fused quartz (ñqz = 1.95 − 0.0 048i), 31 assuming that the sample is being measured at θi = 40 ○ and 65 ○ , respectively.rp and rs are calculated by Eqs. (2) and (3), respectively.Ẽi(ω) is simulated by a normalized double-Gaussian filter with −6 dB cutoff 1. Estimateñ sample (ω) and calculaterp(ω) andrs(ω). 2. Simulate or measureẼi(ω). 3. Add random noise toẼrp andẼrs:  Figure 4(a) shows examples of the simulated detected spectrum with |r| = 1 and |r| = 0.1, respectively, showing a reduced SNR and useful bandwidth with a lower reflectivity. For quartz measurement, 200 groups ofẼrp andẼrs signals are calculated with random noise added, corresponding to 200 groups ofρnoise for each incident angle. Figure 4(b) shows the SNR of tan Ψ at the two incident angles, calculated by mean(tan Ψ)/std(tan Ψ) (mean and std represent average and standard deviation, respectively). Since the value of θB for fused quartz is 62.9 ○ , the SNR of the orange curve calculated at θi = 65 ○ is about one-fourth that of the blue curve, due to the weak p reflection. Interestingly, Figs. 4(c) and 4(d) show that the final errors in n and κ of θi = 65 ○ are less than half the errors obtained for θi = 40 ○ . This verifies our previous discussion that the best sensitivity is achieved around θB and that it has a huge impact on the characterization of our results. At angles far away from θB, small measurement errors can be magnified to the characterization results. The example shows the importance of sensitivity and error analysis in ellipsometry.

B. Angular error
Ellipsometry has a higher angular sensitivity than traditional reflection measurements, mainly due to the need for a large incident angle around θB. This is especially obvious for samples with a large refractive index because, as can be observed in Fig. 2(a), a rapid change in tan Ψ is observed at large angles. Similar to the measurement error, we can estimate the influence of the incident angle error by performing theoretical analysis. Here, we choose bulk high-resistivity silicon (HR-Si), fused quartz, and high- 1.95 − 0.0048i (62.9 ○ ), and 1.54 − 0.01i (57 ○ ), respectively. 31, 32 We assume they are measured under θi = 70 ○ , 65 ○ , and 60 ○ , respectively. Theoreticalρ for the three samples can be calculated by Eqs. (1)-(3). Assuming the actual incident angle is measured as θi + Δθ, the corresponding refractive index obtained from the characterization will have an error Δn compared to the theoretical value, which is shown as a function of Δθ in Fig. 5. A higher angular sensitivity is found for materials with a larger refractive index as expected, but overall, all the three samples are sensitive to the incident angle compared to traditional transmission or reflection measurements. In ellipsometry, attention should be given to the optical set up to reduce Δθ.
The angular accuracy can be improved in two ways. First, good alignment is essential: the THz beams should be aligned parallel to the optical rails. This is because the incident angle is physically measured by the angle between the two rails (or other mechanical components), it is fundamental to ensure that the angle measured represents the angle of incidence. When on-axis lenses are used, properly extending the beam paths on the rails can reduce the error due to the tilting angles of the antennas. Second, the issue of beam divergence should be carefully considered. This issue has been noticed for infrared ellipsometry due to the longer wavelength involved, 1 and it becomes more significant for THz waves. Collimated THz beams provide a near-zero angular spread, however, with a centimeter-level beam diameter as it must be adequately larger than the wavelength to meet the plane-wave approximation. In ellipsometry, the use of large incident angle further enlarges the illuminating area due to elliptical projection, making collimated beams impractical for most samples. This is limited not just by the sample size but also by the surface evenness or filmthickness homogeneity. Focusing the beam introduces a trade-off issue between the spot size and the beam divergence. Physically, the radius (i.e., beam width) of a Gaussian beam w is defined as the distance between the optical axis and the position at which the light intensity drops to e −2 times the on-axis intensity, as shown in Fig. 6. w is related to the beam propagating distance z (in air) as follows: 33 where w 0 is the beam waist, defined as the minimal beam radius achieved at the focal point (z = 0) and λ is the wavelength. As the beam propagates away from the waist, w increases to form an angular spread. Having a smaller w 0 generates a faster variation of w with z, hence, a larger beam spread. This can be expressed mathematically by Eq. (11), which is a hyperbola with an asymptote slope of θ = arctan( λ πw 0 ), as indicated by the blue line in Fig. 6. The expression for θ clearly shows the trade-off issue between the beam divergence and the spot size. Optics with a smaller numerical aperture (NA = sin θ), or equivalently a larger f -number (N = EFL/D, where EFL is the effective focal length, D is the effective aperture), reduce the divergence but increase the focal spot size. To quantitatively show the relationship between the spot size s = 2w 0 and θ, Table II is given as a reference. In the calculation, we assume the collimated THz beam has a diameter (i.e., D) of 30 mm for all frequencies, and it is focused by lenses or mirrors with different EFL values. This assumption gives the condition of w(EFL) = 15 mm. Substituting this into Eq. (11), we obtain w 0 , as well as s and θ for different frequencies. In practice, the sample is placed at the focal point such that only a limited range of the beam will interact with the sample, as illustrated in Fig. 6. Since the beam is less diverging around the focal point, the actual angular spread ϕ is smaller than θ. ϕ can be evaluated by the edge of the beam interacting with the sample, and it is identical to the slope of the tangent at this point, expressed as ϕ = arctan w ′ |z=z 1 . z 1 is the z position at which the edge of the beam interacts with the sample, which is related to the incident angle by tan θi = z 1 w(w 1 ) . Note that the two edges of the beam interacting with the sample have the same spreading direction, thus the total spreading angle is considered as ϕ rather than 2ϕ. Here, we assume θi = 60 ○ to calculate ϕ in Table II. With the defined incident angle, the projecting length of the illuminating area on the sample in the incident-plane direction l is also estimated by l = 2z 1 sin θ i , as shown in Fig. 6. s and l together indicate the minimum sample size required under the specific optics and wavelength. One obvious characteristic observed from the table is that both θ and ϕ are nearly frequency-independent, while s and l are basically proportional to the wavelength. This explains the severity of angular spread issue for longer wavelengths. Using optics with a large EFL offers a smaller beam divergence for all frequencies, while the high-frequency components can be focused to a smaller size. Second, it is found that although both θ and ϕ decrease with increase in EFL, θ is roughly inversely proportional to EFL while ϕ is roughly inversely proportional to EFL 2 . The reason that ϕ decreases much faster comes from the longer depth of focus when using optics with a larger EFL, which means the weakly diverging region near the focus is longer. Note that the effect of beam divergence on the reflection is rather complicated, and one cannot simply equate ϕ and Δθi in Fig. 5 to evaluate or calibrate the error. The divergence deviates the actual reflection away from the plane-wave approximation made in the derivation of Fresnel coefficients. However, the high angular sensitivity indicates the importance of using optics with a smaller beam divergence to reduce the error on the characterization results. We also notice that both s and l are about linearly proportional to EFL. Increasing EFL is more efficient in reducing the divergence ϕ compared to the expansion of the illuminating area. Therefore, a general guideline for the selection of optics is to maximize the EFL as long as the beam spot size can be supported by the sample. Another strategy is to give up some low-frequency components whose illuminating areas are too large for a specific sample.

C. Limited extinction ratio
THz polarizers are mostly wire-grid polarizers (WGPs) made by using subwavelength metallic grids. 34 They can be commercially purchased or self-fabricated by a couple of techniques. In general, they fall into three categories, i.e., bulk-substrate, thin-film, and free-standing polarizers, which have different advantages and limitations. In ellipsometry, the most important parameter is the extinction ratio (ER), defined as the ratio between the transmitted intensities from the passing and blocking directions of a polarizer, expressed as follows: ER = T pass T block . We will later analyze the influence of ER on the measurements. In conventional designs made of single-layer subwavelength metallic grids, ER increases with decreasing wire width and period. The achievable ER is mainly limited by the fabrication technology. ER up to 10 4 over 0.1-2.0 THz with a 140 nm period has been demonstrated. 35 The use of multiple layers of metallic grids is another efficient approach to improve the ER, which normally doubles the ER (in dB scale) compared to the single-layer design. 36,37 Other strategies may also be applied, such as using triangular or sinusoidal metallic surfaces. These methods can be adapted with bulk-substrate or thin-film polarizers. Free-standing polarizers are made using metallic wires such that the ER can only be improved by reducing the wire width and period. 38 However, limited by the mechanical strength, the wire width is typically limited to be above 5 μm. 39,40 As a result, the ER is often insufficient especially in systems with an ultrabroad bandwidth because the ER is basically inversely proportional to the frequency. Nevertheless, free-standing polarizers have two additional attractive characteristics: They provide the best transparency in the passing direction over the three designs as no absorption materials are involved and they are also devoid of Fabry-Perot (FP) effects arising from the multiple reflections within the substrate or the thin film. The second feature is especially useful in time-domain spectroscopy as the multiple reflections either broaden the pulse width to decrease the temporal resolution or introduce pulse echoes in the time domain that can be interfered with reflections from the sample.
In the ellipsometer configuration shown in Fig. 3, the ER of P2 has the largest impact on the measurement accuracy. To show this, we define the ratio between the transmission coefficients of the passing and blocking directions of P2 as follows: It can be seen that ER = |R 2 |. When the ER is not sufficiently high, R cannot be approximated as infinite in practical measurements. In this case, the measured signal when P2 is set in the p direction becomes where the 1/ √ 2 term comes from the projection from the p or s direction of P2 onto the 45 ○ direction of P3, as shown in Fig. 3(b). Ellipsometers usually measure a sample around θB, where |Ẽrp| ≪ |Ẽrs|; hence, the second term in the brackets cannot not be omitted when R is not sufficiently large. For example, if |Ẽrs| = 10|Ẽrp|, |R| > 1000 (i.e., ER > 60 dB) is needed to have the amplitude error forẼ meas rp to be less than 1%. Such a requirement is difficult to achieve for single-layer WGPs, especially at high frequencies where ER decreases. 36,41,42 Figure 7 shows the amplitude and phase of 1/R of a free-standing WGP with a wire width of 5 μm and a period of 12 μm, which has nearly the highest ER among commercially available free-standing WGPs. |R| is less than 100 for frequencies above 0.6 THz. Therefore, in many cases, it is recommended to numerically calibrate the measured signals for the limited ER of P2. To do this, R of P2 should be first characterized as follows: whereẼpass andẼ block are the transmitted electric fields when the passing and blocking directions of P2 are aligned parallel to P1 (and P3), respectively. We rewrite Eq. (13) and add the relationship for the measured s signal as follows: Combining Eqs. (15) and (16) which provide the equations required to calibrate the measured signals using R. Notice that tpass is not measured, but it can take arbitrary values as it will be canceled out when we finally calculatẽ Erp/Ẽrs.The approximation by omitting terms with R 2 is mostly satisfied as R 2 > 100 can be achieved in most polarizers, such that the introduced error is less than 1%. P1 and P3 have the passing direction parallel to the major polarizing direction of the antennas. Assuming that the emitted signal (or the detection sensitivity) is mainly linear, with the emitted electric field (or detector response) in the major polarizing direction asẼmajor and that in the perpendicular direction asẼminor, the emitted field from P1 (or the detected field with P3) becomes E filter = tpass(Ẽmajor +Ẽminor/R). Most antennas have |Ẽmajor| > 10 |Ẽminor|, the influence of the second term will be smaller than 1% with R > 10. Therefore, the ER issue for P1 and P3 is usually ignorable.
Another factor that produces error similar to the ER effect is the orientation error of P2. Analogously, a small angular error for P2 when measuring the p component could lead to a relatively large s component projected onto the passing direction. To set an accurate orientation, all polarizers should be first carefully aligned parallel relative to each other in a transmission configuration, then all of them should be aligned to the sample p − s coordinate in a reflection configuration (e.g., calibrated by a metal mirror), generally referred to as the gravity vertical. In addition, the rotational error should be especially considered when P2 is manually rotated due to the limited the controlling accuracy. It is recommended to reduce the rotating error using a motorized rotational stage such that <0.1 ○ repeatability can be achieved.

D. Pulse shift
Reflection characterization techniques are always sensitive to the phase, including ellipsometry. This is because the maximum phase variation induced by a sample is ±180 ○ . 43 In contrast, the pulse can be delayed for a long distance in transmission to create a huge phase change. 44 Although ellipsometry is a self-reference technique without a phase-uncertainty issue from the sample displacement, other factors such as delay-line positioning error or pulse shift mostly caused by temperature variation in fiber-based THz-TDS systems can generate phase error. Pulse shifts can have a significant influence at high frequencies. The phase error Δϕ is related to the pulse shift τ by Δϕ(ω) = ωτ. For example, a 10 fs pulse shift results in 3.6 ○ , 7.2 ○ , and 14.4 ○ phase errors at 1.0, 2.0, and 4.0 THz, respectively, which can cause significant error in the characterization results at high frequencies. Such an error is very system-dependent, and it is recommended to test both the short-term and long-term (over a few minutes) phase stability by continuously recording a reference signal. Figure 8 plots the measured pulse shifts of two fiber-based THz-TDS systems (determined from the slope of the unwrapped phase), showing very different levels of pulse stability. The variation in system A is found to be correlated with minor variations in the environmental temperature. Here, we introduce a method to calibrate the unknown pulse shift in ellipsometry measurement. 28 The basic principle involves measuring a third signalẼ meas β by rotating P2 to β direction (β ≠ 0 ○ or 90 ○ ). If there is no pulse shift error,Ẽ meas β can be reconstructed fromẼrp andẼrs, mathematically expressed as their linear combination. In contrast, the expression is impossible if there are pulse shifts between the measured p, s, and β signals. We assume there are two pulse shift errors, i.e., τp and τ β , relative toẼrs. This principle allows us to find τp and τ β by minimizing the error betweenẼ meas β and the linear combination of Erp andẼrs. Notice that the ER calibration is not independent of the pulse-shift calibration, and they should be performed together. Table III gives the calibration steps for the pulse shift and ER errors simultaneously.

TUTORIAL scitation.org/journal/app
Step 1 empirically assigns τp and τ β within a searching range according to the potential error of the system. Step 2 calibrates the ER using Eqs. (17) and (18) by taking τp into consideration. Here,Ẽ meas rp andẼ meas rs are again multiplied by √ 2/tpass, but tpass can still be arbitrary as it will be canceled out in step 4.
Step 4 simulates the detected β signal by considering the limited ER effect. tpass in the equation will be canceled out with the scaling factor of 1/tpass contained withinẼ proj β andẼ proj β ′ .
Step 5 compares the simulated and measured β signals by summing up the values in the spectral range with a good SNR, specified by ω 0 and ω 1 , with the potential pulse shift τ β considered. This equation serves as the evaluation function of the algorithm. Theoretically, ΔE β is zero only if both τp and τ β are correctly found. Therefore, τp and τ β are determined at the minimum of ΔE β . Finally, τp found from the minimum can be substituted into the equations in step 2 to calibrate both the ER and pulse shift errors.

E. Other errors
Some other errors may require individual considerations for specific measurements. For example, measuring samples with a wavelength-comparable thickness results in multiple reflections that overlap temporally. At large incident angles, the high-order reflections could have an obvious offset in space, that is, the reflections become only partially overlapped with the surface reflection. The degree of alignment with the detector, and hence the detection sensitivity, is changed compared to the surface reflection. Far sub-wavelength thin films cause negligible changes to the alignment while multiple reflections from bulk substrate can be temporally chopped. Therefore, it is wise to avoid measuring samples with a thickness of around 50-300 μm (depending on the pulse width and sample refractive index). Another potential error is the flatness. Wavelength-comparable roughness, or uneven surface over the illuminating area, scatters or depolarizes the light. Selecting proper optics to control the depth of focus and beam size, or numerically filtering the unaffected spectral range, may reduce these errors.

A. Analytical solution
Characterization refers to extracting the sample optical properties from the measurement, which is the fundamental purpose of spectroscopic ellipsometry. Characterization could be either simple or complicated depending on the sample and the measurement. The simplest case is the analytical solution that is possible when measuring homogeneous and isotropic bulk samples. In this case, combining Eqs. (1)-(3), we have the following: 1 whereεi andεt are the complex permittivities of the incident (i.e., air) and transmitted (i.e., sample) media, respectively, relating toñ byε =ñ 2 . Withεi = 1, Eq. (19) has an analytical solution as follows: n and κ can be determined from εt. For conductive materials, it is convenient to represent the properties by the complex optical conductivityσ, which is related to the permittivity as follows: where ε 0 and ε∞ are the vacuum permittivity and the permittivity, respectively, at high frequencies.

B. Numerical solution and model fitting
Numerical solution usually applies to the cases where only one set of frequency-dependentñ(ω) [or equivalentlyε(ω) or σ(ω)] is unknown but analytical solution is impossible or difficult to be derived. For example, measuring a thin-film sample with the thickness already known, we have the relationship of ρ(n, κ) = tan Ψ exp(iΔ). At each individual frequency, we are mapping (n, κ) to (Ψ, Δ). Therefore, (n, κ) can be found by numerically calculating (Ψ, Δ) over a searching space to find the point that has the smallest difference to (Ψmeas, Δmeas). As the solution can be found independently for each frequency, only two parameters need to be determined at once, thus the computational complexity is low. Many algorithms are available, such as the Newton-Raphson and Gauss-Newton methods that have been used for thin-film characterizations in transmission. 45,46 APL Photonics TUTORIAL scitation.org/journal/app A more complicated situation is when multiple measurements are taken. Measurements at various incident angles, at different sample orientations, or by adding an additional transmission measurement are necessary when the sample has more unknown parameters, such as thin-film samples with an unknown thickness, anisotropic medium, or multiple-layer structures. Actually, even simple isotropic materials for which an analytical or a numerical solution is available benefit from multiple measurements by the improved fitting robustness. 14 This is a common strategy in IR-UV ellipsometry. Redundant measurements are preferred so as to improve convergence and eliminate mathematical ambiguity. However, when there are more unknown variables [such as two sets of n(ω)], the solution becomes a high-dimensional model fitting problem, which significantly increases the computational complexity. Although theoretically, every individual frequency is independent and it is still possible to fit point-by-point to simplify the calculation, it is not recommended as the optimization is not robust when the searching dimension is large while the number of unknown values is not significantly less than that of the measured values. For example, finding [n 1 , κ 1 , n 2 , κ 2 ] by fitting to [Ψ 1 , Δ 1 , Ψ 2 , Δ 2 ] at each frequency is usually not robust and could be trapped into local minimums, leading to unreasonable fluctuations in the optical spectral properties. A more robust and widely adopted protocol is to describe the sample properties by a proper dielectric model and fit to all measured spectra at once. The model essentially reduces the number of unknown values to a few model parameters. For example, let us express a conductive sample by the Drude model (replace i with −i in this equation for the definition ofñ = n + iκ) as follows: 47σ where N, q, m * , and τ are the concentration, charge, effective mass of carriers, and the scattering time, respectively. The permittivity can be further calculated by substituting Eq. (22) into Eq. (21). In the Drude model, typically, only the carrier concentration N and the scattering time τ are unknown. The number of unknown values for the whole spectrum is reduced to two. Using proper dielectric models can significantly improve the fitting convergence to ensure approaching the global minimum. Other commonly used dielectric models include the Lorentz model, Cauchy model, Sellmeier model, and combinations thereof. For complex mixtures, effective medium theories can be used to express the properties as a mixture of two or more inclusions, usually with the volume fractions being unknown, to be found from the fit. These examples indicate another difficulty of performing model fitting-that the selection of the model is highly experience-based, sometimes with trial and error on different models to find the best fit. Recently, machine learning has been proposed for solving this issue by providing an automatic and unambiguous fit to the IR-UV ellipsometric data, 48 which is promising but less meaningful for THz ellipsometry before it is widely available as a standardized characterization technique. Another solution is using a general mathematical model to describe the dielectric behavior. For example, we have formerly shown that the empirical exponential function can provide a good description of n (or κ) for various materials in the THz range as follows: 49 where a, b, and c are three parameters representing the frequencyprofiles of n or κ, hence, leading to six parameters in total to describe the material's complex spectral properties. The purpose of using mathematical models is to reduce the number of unknown variables and characterize the sample unambiguously. Its adaptability is usually limited to materials without resonant features. For broadband ellipsometry, the exponential function may not provide a good approximation for the full bandwidth and more terms may be added to increase the fitness. This method has even been successfully applied to perform ellipsometric characterization of living skin, 11 as shown in Fig. 10. In this work, a double-prism design is adopted to enable two alternative incident angles to be easily switched by moving the prisms vertically. Equation (23) is used to express the complex refractive indices of the stratum corneum (anisotropic) and epidermis.

VI. EXPERIMENTAL DEMONSTRATION
We demonstrate the experimental ellipsometric measurement and data processing of a HR Si wafer (sample A), a low n-doping Si wafer (sample B), and a high n-doping Si wafer (sample C) as examples. The system configuration is the same as that presented in Fig. 3 Table III. The calibration of sample C is shown in Fig. 11 as an example. Figure 11(a) shows the normalized ΔE β calculated by using step 5 in Table III. A single minimum is found at [τp, fs. At this point, the magnitude spectra of theẼ β−P3 ,Ẽ meas β and their difference ΔE β are shown in Fig. 11(b). ΔE β has reached the noise floor, indicating that we have found the minimum available under the system SNR. tan Ψ and Δ with and without calibration are shown in Figs. 11(c) and 11(d), respectively. As the system we used has a high pulse stability, the major difference between them is provided by the ER calibration that corrects the high-frequency signal leakage. However, the singular minimum in Fig. 11(a). and the rapid variation against τp and τ β show the high sensitivity of the algorithm in pulse-shift calibration, which could be very important for systems with low pulse stability.
The characterization was then performed based on the calibrated tan Ψ and Δ. As the samples are bulk and isotropic, Eq. (20) can be applied to obtain the analytical solution. Alternatively, we can apply the Drude model to describe doped semiconductors in the THz range. Therefore, Eqs. (21) and (22) were used to find the Drude parameters by fitting to the whole spectrum ofρ at once for samples B and C. In the model fitting, q = e (electron charge for n-type), ε∞ = 11.7, and m * = 0.26me are used, 50,51 where me is the electron mass. Sample A (thickness 987 μm) and sample B (thickness 538 μm) are sufficiently transparent in the effective bandwidth, hence they are further measured by transmission to evaluate the experimental accuracy. Figure 12 shows the correspondingñ characterized by the different approaches mentioned above for the three samples. A high degree of agreement can be found for the results obtained by different characterization methods, confirming the high measurement and calibration accuracy. Some negative values of κ are found in the analytical solution for samples A and B, which are caused by small phase errors making Δ slightly greater than 180 ○ . Table IV gives the Drude parameters N and τ for samples B and C found from the best fit in 0.2-3.5 THz, as well as the DC (direct-current) resistivities from the model R Drude 0 = 1/σ(0). R Drude 0 matches well with the DC resistivities R 0 measured by the four-probe method. The resistivity can be more accurately determined if the bandwidth can be extended to lower frequencies. To further evaluate the model, we adopt values from the mobility-carrier density curve by Jeon and Grischkowsky. 50 The scattering times from the literature (calculated from the mobility) at the corresponding carrier densities of samples B and C are provided in the table and they highly coincide with the values from the fit. As transmission measurements are limited to carrier densities below 10 17 cm −3 for n-doped Si, 50 sample C is completely opaque to the studied wavelength. In this example, we demonstrate one of the advantages of ellipsometry in accurately characterizing  highly absorptive solids in the THz regime. Similar methods can be applied in other measurements and further combined with multiple configurations (e.g., prism-coupled, different incident angles, sample azimuthal angles), serving as a versatile technique for various sample types.

VII. SUMMARY
THz-TDS has been widely applied as a powerful characterization tool for physical, chemical, and biomedical applications. However, the commonly used transmission and window/prismsupported reflection configurations have limitations in the characterization of thin films, absorptive solids, and complicated structures. THz spectroscopic ellipsometry has a great potential to fill these gaps, despite having a higher requirement in regard to optical alignment, polarization manipulation, and data processing. These technical details, particularly the sensitivity, angular divergence, limited ER, and pulse shift, significantly affect the characterization accuracy, but they are not well known since they are less important in conventional configurations.
In this Tutorial, we quantitatively analyzed errors that commonly occur in THz ellipsometry and their propagation, discussed the reasons, and provided methods for their calibration. Characterization methods are also introduced for different types of measurements. Compared to the conventional ellipsometers at higher frequencies, the major difficulties for THz ellipsometry come from the angular control due to the trade-off between the beam size and the angular spread, as well as the potentially higher phase error in fiber-coupled antennas. Nevertheless, THz ellipsometry also has a few unique advantages. The most promising feature is the coherent generation and detection of THz waves, which significantly simplifies the polarization control. The timedomain sampling scheme also provides a good temporal resolution to easily separate multiple reflections for thick samples, avoiding the out-of-focus and standing-wave issues and simplifying the data analysis. By properly eliminating or reducing the errors, THz spectroscopic ellipsometry can be used as an accurate characterization tool with a very broad adaptability. With the ability to precisely control the incident angle, polarization, and phase, it can also serve as a versatile platform for testing functional devices or for being applied in various interdisciplinary applications.