Gouy shift correction for highly accurate refractive index retrieval in time-domain terahertz spectroscopy

Terahertz spectroscopic measurements are usually performed in focused beam geometry while the standard routine for the retrieval of the sample refractive index assumes plane-wave approximation. In this paper we propose a model for the transmission function which accounts for spatially limited Gaussian terahertz beams. We demonstrate experimentally its validity and applicability for an accurate extraction of the refractive index from experimental data. ©2010 Optical Society of America OCIS codes: (300.6495) Spectroscopy, terahertz; (120. 7000) Transmission. References and links 1. D. Grischkowsky, S. Keiding, M. van Exter, and Ch. Fattinger, “Far-infrared time-domain spectroscopy with terahertz beams of dielectrics and semiconductors,” J. Opt. Soc. Am. B 7(10), 2006 (1990). 2. M. C. Nuss, K. W. Goossen, P. M. Mankiewich, and M. L. O'Malley, “Terahertz surface impedance of thin YBa2Cu3O7 superconducting films,” Appl. Phys. Lett. 58(22), 2561 (1991). 3. M. Misra, K. Kotani, I. Kawayama, H. Murakami, and M. Tonouchi, “Observation of TO1 soft mode in SrTiO3 films by terahertz time domain spectroscopy,” Appl. Phys. Lett. 87(18), 182909 (2005). 4. C. Kadlec, F. Kadlec, H. Němec, P. Kužel, J. Schubert, and G. Panaitov, “High tunability of the soft mode in strained SrTiO3/DyScO3 multilayers,” J. Phys. Condens. Matter 21(11), 115902 (2009). 5. 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). 6. L. Duvillaret, F. Garet, and J.-L. Coutaz, “A Reliable Method for Extraction of Material Parameters in Terahertz Time-Domain Spectroscopy,” IEEE J. Sel. Top. Quantum Electron. 2(3), 739–746 (1996). 7. I. Pupeza, R. Wilk, and M. Koch, “Highly accurate optical material parameter determination with THz timedomain spectroscopy,” Opt. Express 15(7), 4335–4350 (2007). 8. S. Feng, H. G. Winful, and R. W. Hellwarth, “Gouy shift and temporal reshaping of focused single-cycle electromagnetic pulses,” Opt. Lett. 23(5), 385–387 (1998). 9. P. Kužel, M. A. Khazan, and J. Kroupa, “Spatio-temporal transformations of ultrashort terahertz pulses,” J. Opt. Soc. Am. B 16(10), 1795–1800 (1999). 10. M. T. Reiten, S. A. Harmon, and R. A. Cheville, “Terahertz beam propagation measured through threedimensional amplitude profile determination,” J. Opt. Soc. Am. B 20(10), 2215 (2003). 11. H. Kogelnik, “On the propagation of Gaussian beams of light through lenslike media including those with a loss or gain variation,” Appl. Opt. 4(12), 1562 (1965). 12. A. Dreyhaupt, S. Winnerl, T. Dekorsy, and M. Helm, “High-intensity terahertz radiation from a microstructured large-area photoconductor,” Appl. Phys. Lett. 86(12), 121114 (2005). 13. J. Dai, J. Zhang, W. Zhang, and D. Grischkowsky, “Terahertz time-domain spectroscopy characterization of the far-infrared absorption and index of refraction of high-resistivity, float-zone silicon,” J. Opt. Soc. Am. B 21(7), 1379–1386 (2004). 14. J. Petzelt, P. Kužel, I. Rychetský, A. Pashkin, and T. Ostapchuk, “Dielectric response of soft modes in ferroelectric thin films,” Ferroelectrics 288(1), 169–185 (2003). 15. C. Kadlec, V. Skoromets, F. Kadlec, H. Němec, J. Hlinka, J. Schubert, G. Panaitov, and P. Kužel, “Temperature and electric field tuning of the ferroelectric soft mode in a strained SrTiO3/DyScO3 heterostructure,” Phys. Rev. B 80(17), 174116 (2009).


Introduction
The time-domain terahertz (THz) spectroscopy is now a well established method for measurements of optical constants of bulk [1] and thin-film samples [2][3][4] in equilibrium over a wide range of frequencies extending from a few tens of gigahertz up to a few terahertz.The technique requires measuring the temporal profile of the electric field of a THz pulse (wave form) transmitted through the investigated sample.The complex spectrum of the wave form is normalized by a reference spectrum (obtained when the sample is removed from the THz beam path or when it is replaced by a substrate with known properties); this ratio provides the complex transmission function of the sample.The accuracy of the determination of the complex refractive index in the THz spectral range then essentially depends on the precision of the sample thickness determination [5].
Frequently, the temporal windowing is applied to the time-domain data.In this case the wave form is divided into several intervals, each of them carrying a single pulse belonging to different orders of internal Fabry-Pérot reflections in the sample.The complex refractive index of the sample is then retrieved from these "echoes" [6].A separate analysis of the echoes allows one to find the optimum thickness of the samples from the THz spectra [5], which, in the case of thin films on substrate, may substantially increase the accuracy of the retrieved dielectric spectra of the film [4].Recently an improvement of the signal processing has been reported [7] allowing one to achieve a metrological precision in the THz timedomain spectroscopy of bulk samples.
Nevertheless, it has not been extensively discussed yet in the THz spectroscopic literature, that the accuracy of the THz refractive index depends also on the correctness of the model transfer function serving for the refractive index retrieval.Usually, the plane-wave approximation is used to relate the measured complex transmission function to the complex refractive index of the sample.However, it has been shown [8,9] that spatially limited ultrashort THz pulses change their phase dramatically if they are strongly focused; this phase shift is called the Gouy shift.The tight focusing is a frequent issue in the THz spectroscopy as the typical dimensions of the samples exceed the sub-THz wavelengths only by a factor smaller than or comparable to 10.It follows that the next step in the improvement of the accuracy of THz spectroscopy must take into account the spatial properties of the THz beam.
In this paper we go beyond the plane-wave description in analyzing the transmission function of bulk samples and we show that the beam focusing is responsible for a small but clearly measurable systematic deviation of the refractive index.We suggest corrections of the transfer function serving for refractive index retrieval and for the sample thickness determination from the THz data.

Theoretical description
In this paper we describe the THz beam as a Gaussian TEM 00 one.In the literature one finds measurements showing the possible presence of higher-order Gaussian modes in the THz beam [10]; however, later it will become clear that the magnitude of the effect is such that for spectroscopic purposes a simplified description is sufficient.We consider an experimental setup as depicted in Fig. 1, where the emitted THz beam propagates in vacuum and it is focused first onto a sample and subsequently onto a sensor (by a transformation optics with focal length f).We assume that the distances between mirrors or lenses and focal points are carefully aligned when the sample is removed from the beam path.By introducing a homogeneous sample with thickness d and refractive index n, the middle (sample) focus position is shifted by Δ = d(n -1)/n in the z direction (note that in general n may be a frequency-dependent quantity and, consequently, the shift may slightly vary with frequency).The focus position at the sensor is then displaced in the same direction by (see Fig. ( ) within the approximation Δ << |a -f | which is usually well fulfilled.Within the plane-wave approximation, the frequency component ν of the THz field reads: where e 0 denotes the THz wave amplitude and 2 /c k πν = is the wave vector of the radiation.The transmission function of the sample (normalized by a phase factor of propagation in vacuum during the reference measurement) then reads: where N = n + iκ is the complex refractive index of the sample, the integers M 2 ≥ M 1 ≥ 0 define the windowing conditions of the time domain signal and ( ) ( ) Within the Gaussian-beam approximation, the on-axis THz field of the spectral component ν close to the focus (located at z = 0) reads [11]: ; e x p w e z e ikz i z w z where the on-axis amplitude variation is reduced due to the expansion of the beam width out of the focus ( ) and the Gouy shift alters the phase of the field.These expressions make use of the Rayleigh confocal parameter Obviously the on-axis THz field is controlled solely by the beam-waist size w 0 .The amplitude in (5) displays a quadratic dependence on the focus displacement at the sensor and this contribution is expected to be small.At the end we will estimate its magnitude.The largest contribution of the beam focusing to the transmission function comes from the linear phase term defined by the Gouy shift: where ϕ was introduced in (4) and the dimensionless frequency-dependent parameter β is defined as: ( ) Provided that β(ν) is determined for a given experimental setup, the contribution of the Gouy shift to the transmission function can be quantified.The phase change of a focused beam (including the Gouy shift contribution) related to the direct pass through the sample then reads: 1 .
When evaluating the phase of the m-th Fabry-Pérot internal reflection in the sample, we must take into account that folding of the partially reflected beam inside the sample is equivalent to an increase of the distance between the focal plane and the subsequent focusing element by 2md, i.e.
( ) and ( ) The total phase change of the m-th echo is The correction of the amplitude is then deduced from ( 6): ( ) The formula (3) can be used for the evaluation of the refractive index of the measured sample ϕ , and G m A are defined by ( 11), (14), and (15), respectively.

Influence of the Gouy shift on the refractive index determination
In order to estimate the impact of the term β on the evaluation of the refractive index within the plane wave approximation, let us assume for a short while that the absorption of the sample is small: κ << n.In this case, the phase change due to the Fresnel transmission coefficient is very small and can be neglected.Let n be the true refractive index of the material (this value would agree with that found in the collimated-beam geometry) and n 1 be the value which one would find by usual evaluation method in the focused-beam geometry, i.e., by using (2).We assume that the time windowing of the experimental waveform is possible and we consider the direct pass of the beam through the sample.The refractive index is determined by the phase of the wave (5) which, taking into account ( 9), yields: The left-hand side of (17) assumes that the experimental phase change is entirely due to the propagation through the sample (standard evaluation method using the plane wave approximation and leading to n 1 ) while the right-hand side attributes an appropriate part of the phase shift to the Gouy shift.Equation ( 17) predicts that the determination n 1 overestimates the true value n of the refractive index.The correction of the refractive index due to the Gouy shift is relatively small, of the order of β: note that this correction does not depend on the sample thickness.One can easily show that if the Fresnel transmission terms were included in (17), an additional correction in (18) would appear of the order of ~β κ/n 1 which is negligible.The Fresnel transmission terms then can be safely omitted in this analysis of weakly absorbing materials.
Similarly, if we consider the m-th echo of the THz signal, the phase of the wave (5) should be expressed by using (13); the correction of the refractive index calculated from this echo then reads: ) The correction term is thus different for each echo, namely it can change the sign for higher order echoes.It follows that the Gouy shift has also an influence on the sample thickness determination from several echoes of time-windowed signal.This method has been first described in [5] and relies on the fact that the thickness of the medium probed by the m-th echo varies with the echo number m: (2m + 1)d.The retrieved refractive index then depends on the echo number if a wrong thickness of the sample is used in the retrieval procedure.Our findings show that this method needs to be revisited, too, if a high precision of the thickness determination is required.Therefore the parameter β appears crucial for the evaluation of the Gouy shift correction and its magnitude determines directly the importance of the correction.

Determination of parameter β
The parameter β should be determined experimentally for each particular THz spectroscopic setup.It would be possible, using (10) and ( 8), to determine β by measuring the beam waist size for each frequency component.However, it is much more convenient and precise to evaluate β from the measured Gouy shift.In fact, we can estimate how much the focus is shifted near the sensor due to the inserted sample and we may compensate for this displacement approximately by shifting the sensor in the same direction by: where n (≈n 1 ) is an average THz refractive index of the sample.
In order to determine the parameter β we suggest an experiment made of 3 consecutive measurements.(1) Standard measurement of a reference wave form without the sample; (2) inserting the sample and measuring the sample wave form: determination of n 1 within the plane wave approximation (the beam focus is displaced by Δ' with respect to the sensor); (3) second measurement of the sample wave form, now with the sensor displaced by 2D (i.e. the beam focus is shifted by -Δ' with respect to the sensor), determination of n 2 within the plane wave approximation: By comparing (17) and (21) we find: ) This formula can be used for the determination of β with a good precision and within the linear approximation (small Δ'/z 0 ).Also, β(ν) should not depend on the properties of the particular sample.Repeated measurements with several samples (various thicknesses and refractive indices) should yield the same β; this constitutes a verification of the validity of our model.

Synopsis
In order to verify the Gaussian-beam model described in the previous sections several kinds of experiments may be envisaged.At the first sight it might seem useful to characterize a sample in transmission by using a large collimated THz beam and, in the second step, to add two additional lenses into the setup and make an additional measurement with the sample placed into a tight focus.However, we deal here with very tiny effects, so that this experiment would be prone to errors stemming from the following two issues.(1) A large sample (at least 6 cm in diameter) is needed for such experiment to be reliable.It is then necessary to ensure a perfect homogeneity in its dielectric properties and, above all, in its thickness over the whole section.
(2) Inserting the optics into the setup is required between the two connected experiments.This step may lead to misalignment of the beam axis.Anyway the function β(ν) would have to be experimentally determined for the two setups.
For these reasons we find more useful to characterize samples in two spectrometer configurations where either a moderately or a more tightly focused beam is used.The fact that we use two different states of the focused beam gives us the possibility to evaluate the contribution of the Gouy shift.Comparable areas of the samples are probed in this case and, moreover, the variation of the THz beam waist is achieved only by a different focusing of the optical beam onto the emitter without any change in the layout of the THz optics.We proceed in our verification in three steps: 1) The function β(ν) is determined for a large number of samples for the two configurations.The fact that β(ν) is found to be sample independent constitutes the first indirect proof of the model.This part of the study provides us with all the data required for the Gouy shift evaluation.
2) We characterize a single sample within the two setups and we show that the retrieved refractive indices differ considerably if the Gouy shift is neglected while the agreement is significantly improved when it is taken into account.
3) In section 3 we have shown that the Gouy shift contribution differs for different echoes (Fabry-Pérot internal reflections in the sample).The most precise verification of the Gouy shift contribution then consists in a measurement of a direct pass and the first echo transmitted through a thick sample.In this case two spectra of the refractive index are obtained from a single measurement without any intermediate adjustment of the experimental setup.The quality of the match of these spectra attests the model.

Experiments
Our experimental geometry is shown in Fig. 2; the THz pulse transmitted by an emitter is focused onto the sample by an ellipsoidal mirror and another identical ellipsoidal mirror focuses the beam onto the sensor ([110]-oriented 1 mm thick ZnTe electro-optic crystal).The focusing of the beam within this setup can be schematically depicted as in Fig. 1 where the focal length of both mirrors is f = 7.5 cm and a = a' = 2f.In this case (a/f -1) 2 = 1 in the above equations, namely in Eqs.(10) and ( 22) which specify the parameter β.The experimental setup can be operated with two different laser sources and emitters, respectively.The first configuration (denoted as Setup I) uses a femtosecond oscillator (Coherent, Mira); its pulse train is moderately focused (~1.2 mm diameter) onto an interdigitated photoswitch TeraSED [12] which emits broadband THz pulses (0.1-3 THz).The second one (Setup II) employs a femtosecond amplifier (Quantronix, Odin); its collimated beam excites a 1 mm thick [110] ZnTe crystal where it generates the THz pulse (0.15-2.5 THz) by optical rectification.The initial diameter of these two THz beams is different and this leads to different β(ν) for the two configurations.This allows us to measure samples with a more or less tightly focused THz beam, compare the results and conclude on the validity of our model.For practical reasons we performed displacements of the emitter in the third step of the experiment proposed in the previous section; these are equivalent to the same displacements of the sensor in our setup.We selected a number of samples with various values of average refractive index in the THz range and with various values of their thickness; their properties are summarized in Table 1.
The experiments were performed as described above; in Fig. 3 we show the measured spectra of the parameter β for both setups.One observes that for each of the two setups, the individual curves β(ν) determined with various samples coincide within experimental noise.This result validates our theoretical analysis.From the magnitude of β we estimate that the expected correction of the refractive index is ~0.01 in the low-frequency part of the spectrum and it practically vanishes close to 2 THz.The correction of the amplitude determined by the coefficient A m does not exceed 0.25% for any of the investigated samples for m = 0, 2 and it is even much smaller for m = 1.This is smaller than the typically achieved accuracy of the amplitude measurement in time-domain THz experiments, which means that the amplitude correction can be neglected in a majority of cases encountered in practice.The correction may become significant only for very thick samples (e.g. for a 6 mm thick Si wafer one obtains A 0 = 0.976 in Setup I at 0.25 THz).It is possible to calculate the beam-waist radius w 0 and the confocal parameter z 0 from an average β for both setups by using ( 10) and ( 8) (see insets in Fig. 3).These results show that in Setup I we work with more tightly focused THz beam which presents an advantage in measuring smaller samples; however, the Gouy shift correction is higher in this case.In [9] we show that the long-wavelength part of the broadband THz pulses may be partially cut by the mirror clear aperture A due to the beam diffraction.This phenomenon occurs namely for frequencies below the cut-off frequency ν m , which depends on the geometrical arrangement [9]: where L is the emitter-mirror distance and w E is the beam-waist size at the emitter.The values of the cut-off frequency for our setups are shown by dotted lines in the insets of Fig. 3 and explain an increase of w 0 and z 0 at the low-frequency part of the spectrum.

Discussion
We address here some aspects of the refractive index determination of the thickest samples of the investigated series; they exhibit the largest Gouy shift correction.The BaF 2 sample was investigated in both experimental setups; it is then possible to compare the refractive indices retrieved from these two experiments within the two discussed models (plane-wave and Gaussian approximation).In Fig. 4 we show the low-frequency part of the refractive index of this material retrieved within the plane-wave approximation and the Gaussian beam approximation from experimental data obtained in Setups I and II.All four curves converge towards the same behavior at higher frequencies.In setup I the beam is more tightly focused and therefore the plane-wave approximation provides a clearly higher refractive index than that calculated within the frame of Gaussian beam optics; the difference between these two curves does not vanish even above 1 THz.Comparison of the "true" refractive indices obtained for the two setups, i.e. those calculated within the Gaussian beam approximation, shows that they match nicely above 0.5 THz and we observe a small difference (< 0.002) below this frequency.The discrepancy is probably related to the accuracy of the determination of β at these low frequencies where the data in Fig. 3 show quite a noisy behavior.We performed a detailed measurement of the refractive index of a thick wafer of Si in Setup I. Data corresponding to the direct pass and to the first echo in the time domain were acquired and processed within both the plane-wave and Gaussian beam approximations.At the same time we performed a wafer thickness optimization by using a procedure similar to that reported in [5] with an appropriate transmission function.The results are shown in Fig. 5.The nominal resistivity of the wafer is ρ = 50 Ω cm so that a weak interaction of the THz pulse with free carriers was observed.Figure 5(a) displays the refractive and absorptive indices retrieved by using the Gaussian beam model.The presence of conduction band carriers is marked by a decrease of the refractive index and increase of the absorptive index in the low-frequency part of the spectrum.The agreement of the data from direct pass and from the first echo is very satisfactory; also the fit of the spectra by the Drude function is excellent and provides the background permittivity of silicon ε ∞ = 11.66 (refractive index 3.415), concentration of carriers n c = 5 × 10 13 cm −3 and momentum scattering time τ ≈350 fs (yielding ρ = 46 Ω cm in agreement with the nominal value); the optimized sample thickness is d = 1.9326 ± 0.0002 mm in agreement with the mechanical measurement (1.933 ± 0.002 mm).
The value of 3.415 of the background refractive index of our Si wafer seems to be in a quite good agreement (difference smaller than 1 per mille) with the so far most precise value of 3.4175 reported for high-resistivity silicon [13].In fact, the authors of [13] cautiously optimized experiments in order to obtain the highest possible precision in the determination of the refractive index of Si, in particular they used a significantly higher thickness than that of our sample which decreases the surface layer/roughness effects and makes their measurements more precise.The influence of surface, expected to play a greater role in our case, may lead to a slight underestimation of the retrieved refractive index if the thickness optimization procedure is used.Although no explicit details about the waist size inside the sample in [13] were provided, it is possible to estimate from the sample dimensions that the beam radius was of the order of 1-2 cm, which provides the Gouy shift correction of the refractive index of approximately 0.5-1 × 10 −4 , a value consistent with the observed decrease in the refractive index in Fig. 4  Figures 5(b,c) show the results of the retrieval within the plane wave approximation.In Fig. 5(b) we used the same thickness as retrieved previously (d = 1.9326 mm) for the calculation of the refractive index.This leads to the same value of ε ∞ for the two data sets (direct pass and first echo); by contrast, the discrepancy at low frequencies is quite high due to the missing Gouy shift contribution in the model.In Fig. 5(c) we proceeded to the thickness optimization procedure within the plane-wave model, which shifts the refractive index to values larger than 3.42.Nevertheless, the discrepancy between the two data sets is still clear.The optimum thickness in this case is 1.926 mm.
We have shown that the Gouy shift correction of the plane-wave transmission function may have significant consequences for an accurate determination of the refractive index and for the determination of the sample thickness from THz time-domain measurements.This may be a very important issue for the measurements of thin films.It is known [14] that the accurate knowledge of the substrate refractive index n s and thickness d s is crucial for the determination of THz dielectric spectra of the film.The error in the refractive index Δn of the thin film in terms of these quantities reads: where d f is the film thickness and d r is the thickness of the substrate serving for the reference measurement.Taking into account the small value of d f (frequently sub-micron) the issue of precise determination of the thickness and refractive index of a bulk substrate becomes critical.Recently, we published papers [4,15] where we pay attention to determine accurately the dielectric response of strained SrTiO 3 films with a total thickness of 100-200 nm; the Gouy shift correction exposed in the current paper has been taken into account in the procedure for substrate optimization.

Conclusion
In this paper we show that in the common geometry used for THz time-domain experiments with a focused THz beam, the Gouy shift occurring along the beam propagation path has a non-negligible influence on the measurement results.This effect has remained, up to now, largely disregarded.While the correction is often small, in some cases, the procedure we describe may bring an important improvement in the precision of the refractive index retrieved from the measured data.Namely, the correction spectral function β -representing the tightness of the THz beam focusing -may be high enough to give rise to clearly noticeable distortions of the spectra, especially in the low-frequency part.This correction is particularly important for the investigation of materials exhibiting a flat refractive index spectrum with a low imaginary part and with a weak feature superimposed at the lowest frequencies (see the above example of weakly doped silicon).In measurements of thin films, an accurate determination of the substrate refractive index spectra taking into account the Gouy shift may improve dramatically the precision of the thin film spectra, notably if there is no reference substrate of the same thickness available.Also, in the case of optically thick bulk samples, the described approach allows one to eliminate a source of systematic error in determining the sample thickness with a high precision, including metrological applications.

Fig. 1 .
Fig. 1.Scheme of an experimental setup containing 3 focal planes (for an emitter, sample and sensor).The displacement Δ of an object plane yields a displacement Δ' of its image.The arrows indicate the shifts of foci.

Fig. 5 .
Fig. 5. Refractive index of a 1.93 mm thick Si wafer calculated from the direct pass (closed symbols) and from the first echo (open symbols).Retrieval by using (a) Gaussian beam approximation and sample thickness optimization (d = 1.9326 mm); (b) Plane wave approximation with a thickness of d = 1.9326 mm; (c) Plane-wave approximation with a thickness optimization (d = 1.926 mm).

Table 1 . Summary of measured samples. The thickness of samples is determined by using a mechanical gauge.
[13]3].
Fig. 4. Refractive index of BaF 2 retrieved from experiments with both setups by using planewave approximation and Gaussian beam approximation.