## Abstract

We develop a generalized causality-based framework for determining the time origin in terahertz emission spectroscopy. Our framework is formulated in terms of a multiply subtractive Kramers-Kronig relation and can treat all major mechanisms of terahertz emission, which include the occurrence of a delta-function-like instantaneous polarization observed typically in nonlinear optical processes. We show that a function derived within our framework properly determines the positions of *t* = 0 both for simulated terahertz waveforms and for a measured one obtained in biased conjugated polymers. This function will be useful for an in-depth understanding of ultrafast phenomena involving terahertz emission in various optoelectronic materials.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Causality is a fundamental concept in physics and is known to regulate the relationship between the real and imaginary parts (or between the amplitude and phase parts) of optical spectra through the Kramers-Kronig (K-K) transformation [1–3]. It has been reported that causality provides a practical method of compensating for a partial lack of spectral information in terahertz (THz) time-domain spectroscopy, where K-K relations can be used together with the values of both amplitude and phase known at some reference frequencies (called anchor points) within the measurement frequency range. Since the singly subtractive K-K relation was successfully exploited by Lucarini *et al*. to remove the artificial phase shift in THz reflection spectroscopy due to a spatial misplacement between the sample and reference materials [4], it has furthermore been applied to THz emission spectroscopy [5] and THz transmission spectroscopy [6,7].

In a previous paper [5], we proposed a causality-based method for determining the time origin (*t* = 0) in THz emission spectroscopy and demonstrated that a misplacement of the time origin can be detected with an accuracy that is an order of magnitude higher than the given temporal resolution of THz electric field. It has been shown that such an accurate time origin in THz emission allows us to obtain an in-depth physical understanding of ultrafast charge transport in semiconductors [8,9]. However, the core theoretical framework of the causality-based method was limited to situations where THz emission does not include the contribution of a delta-function-like instantaneous polarization, which is widely observed in nonlinear optical processes such as optical rectification [10].

In this work, we developed a generalized causality-based framework that can determine the time origin for all major mechanisms of THz emission including the occurrence of a delta-function-like instantaneous polarization. A multiply subtractive K-K relation [11,12] was used to derive a function for evaluating how well the causality in THz emission is satisfied versus a possible time-origin misplacement *δt*. We found that the function reveals the exact positions of *t* = 0 for trial THz waveforms and furthermore works for a recent example of THz waveforms emitted from biased conjugated polymers.

## 2. Generalized formulation of causality in THz emission spectroscopy

We provide the formulation of causality in transient THz emission phenomena and derive the function for determining the time origin, i.e., the moment at which real/virtual electron-hole pairs created by the peak of an ultrashort optical pulse start emitting a THz electromagnetic wave. Let us consider a general situation where the emitted THz electric field *E*(*t*) has conceivable retarded and instantaneous components and can be expressed as

*A*–

*C*are real-valued coefficients. Causality requires that any advanced response be absent and therefore

*E*

_{ret}(

*t*) = 0 for

*t*< 0. Note that Eq. (1) is a generalized version of the expression for

*E*(

*t*) reported previously in [5] and allows us to treat optical rectification, which is known to be one of the major mechanisms of THz emission [10].

By using the theory of complex functions in a way similar to that of [5], we obtain

*Ẽ*(

*ω*) is the complex Fourier transform of

*E*(

*t*) defined by

*A*has disappeared from Eq. (4) owing to a basic property of the Cauchy principal value. The other unknown coefficients

*B*and

*C*can also be eliminated analytically from Eq. (4) when both Re

*Ẽ*(

*ω*) and Im

*Ẽ*(

*ω*) are given at one or more reference angular frequencies (i.e., anchor points), as shown below. We would like to emphasize that the present formulation is efficacious in the elimination of

*C*, which cannot be achieved even if Eq. (1) is adopted in [5]. Setting two anchor points of

*ω*=

*ω*

_{0}and

*ω*

_{2}, we obtain

*E*(

*t*) is a real-valued function and hence Re

*Ẽ*(–

*ω*) = Re

*Ẽ*(

*ω*). More anchor points were introduced into the present formulation than into the previous one [5] so that they could cover various spectral shapes of

*Ẽ*(

*ω*) more properly. Equation (5) is a kind of doubly subtractive K-K relation, which is known to provide faster convergence for the integral over the semi-infinite frequency range than the conventional K-K relation does [11].

Now, we deal with the case where a time-origin misplacement *δt* may accompany the measurement of THz electric field *E*(*t*). Because it gives rise to the artificial shift by –*ωδt* in the measured phase spectrum *θ*(*ω*) = arg*Ẽ*(*ω*) [5,8], Re*Ẽ*(*ω*) and Im*Ẽ*(*ω*) in the formulas derived above should be replaced by *ρ*(*ω*)cos*θ*_{cr}(*ω*) and *ρ*(*ω*)sin*θ*_{cr}(*ω*), respectively, where *ρ*(*ω*) = |*Ẽ*(*ω*)| is the amplitude spectrum and *θ*_{cr}(*ω*) = *θ*(*ω*) + *ωδt* is the corrected phase spectrum. The use of Eq. (5) at another anchor point *ω* = *ω*_{1} leads to the condition for *δt* that *K*(*δt*) = 0, with *K*(*δt*) defined by

*δt*-derivative of the causality-based function

*K*(

*δt*) in Eq. (6) corresponds to the previous version of

*K*(

*δt*) in [5], except for the different numbers of anchor points used in the present and previous formulations. Equation (6) is therefore expected to have useful properties for finding the actual solution of

*δt*numerically. Below,

*ω*

_{0}and

*ω*

_{2}will be set to the quarter maxima of the amplitude spectrum

*ρ*(

*ω*) by following the manner suggested in [5], while

*ω*

_{1}will be set to the maximum or half maxima of

*ρ*(

*ω*).

## 3. Results and discussion

#### 3.1 Verification with simulated THz waveforms

To examine whether the generalized framework described above works properly, we computed the causality-based function *K*(*δt*) for trial THz waveforms simulated with no time-origin misplacement. Figure 1(a) shows a trial THz waveform *E*(*t*) simulated for a transient current of *J*(*t*) = *J*_{0}Θ(*t*)exp(–*γt*)cosΩ*t*, to which the previous version of *K*(*δt*) was applicable [5]. Here, *J*_{0} is the magnitude of current, Θ(*t*) is the unit step function, Ω is the resonance angular frequency, and *γ* is the damping rate; the numerical data on *E*(*t*) was prepared with parameters of Ω/2*π* = 1.5 THz and *γ* = 1.1 THz and was furthermore convolved with a system response function that gives a temporal resolution of *τ*_{res} = 0.30 ps. This type of *E*(*t*) has been observed for Bloch oscillations in semiconductor superlattices [8,13] and gives an example of *A* ≠ 0 and *B* = *C* = 0 in Eq. (1). Figure 1(b) shows the amplitude spectrum *ρ*(*ω*) (solid curve) and the phase spectrum *θ*(*ω*) (dashed curve) calculated by performing the fast Fourier transform of *E*(*t*) in Fig. 1(a); *ρ*(*ω*) has the resonance peak at *ω*/2*π* = 1.5 THz, while *θ*(*ω*) has a change by ~*π* across this frequency.

We fed the numerical spectral data *ρ*(*ω*) and *θ*(*ω*) into Eq. (6) together with anchor points *ω*_{0}, *ω*_{1}, and *ω*_{2}. The causality-based functions *K*(*δt*) computed with (*ω*_{0}/2*π*, *ω*_{1}/2*π*, *ω*_{2}/2*π*) = (0.98, 1.22, 2.06), (0.98, 1.50, 2.06), and (0.98, 1.79, 2.06) THz in an integration frequency range of *ω'*/2*π* = 0.0–4.0 THz are shown in Fig. 1(c) by curves 1–3, respectively. As seen in the figure, the present generalized version of *K*(*δt*) indeed has properties similar to those of the previous version [5]: *K*(*δt*) fluctuates greatly for *δt* < 0 (because causality is violated if the time origin is shifted forward from the current position) but is nearly equal to zero for *δt* > 0 with only a small fluctuation (because causality still holds if the time origin is shifted backward). We find that the small fluctuation, which arises from a finite temporal length of *E*(*t*) [5], is damped the most at *δt* ~0.30 ps and the deviation from the baseline (i.e., the *K* = 0 line) increases significantly as *δt* decreases from there, particularly in the magnified view of curves 1 and 3 in Fig. 1(d) (curve 2 exhibits it less clearly but is consistent with curves 1 and 3). This behavior reflects the fact that the THz signal in Fig. 1(a) penetrates down to *t* ~–0.30 ps owing to the given finite temporal resolution of *τ*_{res} = 0.30 ps. Thus, we can obtain the actual value of time-origin misplacement *δt* by following the same procedure as reported in [5]: we first search for the point *δt'* at which *K*(*δt*) starts deviating from the baseline versus possible misplacement *δt*, and we then estimate the actual misplacement to be *δt* = *δt'* – *τ*_{res}.

Here, we provide a practical criterion for specifying the value of *δt'* more quantitatively. When we consider the largest neighboring peak of *K*(*δt*) on the lower side of *δt* ~0.30 ps, i.e., the negative peak of curve 3 at *δt* = –0.23 ps in Fig. 1(d), we find that curves 1 and 3 at *δt* ~0.30 ps are as small as 1.5% of that peak height in terms of absolute values and we can determine that *K*(*δt*) starts deviating from the baseline at *δt'* = 0.31 ± 0.03 ps. Thus, the actual misplacement is given by *δt* = *δt'* – *τ*_{res} = 0.01 ± 0.03 ps, in which range the correct solution (i.e., *δt* = 0 ps for this simulation) is included properly. The same criterion applies to the other examples shown below. Because the clarity in the behavior of *K*(*δt*) around *δt'* depends on the choice of anchor points, we recommend preparing at least three sets of anchor points to obtain the value of *δt* together with its uncertainty.

Figure 2(a) shows another trial THz waveform *E*(*t*) simulated for a delta-function-like instantaneous polarization with a temporal resolution of *τ*_{res} = 0.030 ps. This type of *E*(*t*) gives an example of the important additional term with *C* ≠ 0 in Eq. (1). As seen in the figure, *E*(*t*) has a symmetric shape with respect to the largest negative peak at *t* = 0. Figure 2(b) shows the amplitude spectrum *ρ*(*ω*) and the phase spectrum *θ*(*ω*), which were obtained from *E*(*t*) and fed into Eq. (6). The causality-based functions *K*(*δt*) computed with (*ω*_{0}/2*π*, *ω*_{1}/2*π*, *ω*_{2}/2*π*) = (5.63, 8.50, 33.9), (5.63, 17.7, 33.9), and (5.63, 28.9, 33.9) THz in an integration frequency range of *ω'*/2*π* = 0.0–50.0 THz are shown in Fig. 2(c) by curves 1–3, respectively. The deviation in *K*(*δt*) from the baseline starts at *δt'* = 0.029 ± 0.003 ps [see also the magnified view in Fig. 2(d)], indicating that the actual misplacement can be estimated properly to be *δt* = *δt'* – *τ*_{res} = –0.001 ± 0.003 ps.

The same analysis as made in Fig. 2 also works for an instantaneous polarization with an order-of-magnitude lower temporal resolution of *τ*_{res} = 0.30 ps treated in Fig. 3. The causality-based functions *K*(*δt*) computed with (*ω*_{0}/2*π*, *ω*_{1}/2*π*, *ω*_{2}/2*π*) = (0.56, 0.85, 3.40), (0.56, 1.77, 3.40), and (0.56, 2.89, 3.40) THz in an integration frequency range of *ω'*/2*π* = 0.0–5.0 THz are shown in Fig. 3(c) by curves 1–3, respectively. The deviation in *K*(*δt*) from the baseline starts at *δt'* = 0.31 ± 0.04 ps [see also the magnified view in Fig. 3(d)], revealing the actual misplacement properly to be *δt* = *δt'* – *τ*_{res} = 0.01 ± 0.04 ps. Putting together the three different examples provided in Figs. 1–3, we find that the uncertainty in determining the time origin from the behavior of *K*(*δt*) is typically ± 10% of the given temporal resolution *τ*_{res}.

Now, let us make a few remarks on the present method. First, it excels the maximum entropy method (MEM) [14] at treating the additional term with *C* ≠ 0 in Eq. (1); the MEM made wrong predictions about the time origin in Figs. 2(a) and 3(a) because the MEM is not good at deriving flat phase spectra [see Figs. 2(b) and 3(b)] from amplitude spectra. Second, the present method is expected to work properly even for the cases where *E*(*t*) in Eq. (1) has several nonzero terms simultaneously; we treated separate cases for simplicity because Eq. (5) consists of linear operations and includes the principle of superposition for *E*(*t*).

#### 3.2 Application to measured THz waveforms

Here, we demonstrate how the present method can be applied to a recent example of actual THz waveforms emitted from biased conjugated polymers. Figure 4(a) shows the THz waveform *E*(*t*) observed for a polythiophene film under an in-plane bias electric field of 12 kV/cm [15], with the time origin set to an arbitrary position. The THz electric field was detected using a 0.2-mm-thick (110) ZnTe electro-optic sensor, whose response function calculated within a theoretical model [16,17] had a nearly flat spectral shape up to 2 THz and thus led to a temporal resolution of *τ*_{res} = 0.50 ps. As seen in the figure, *E*(*t*) has the largest negative peak at *t* = 2.0 ps and smaller positive peaks. The amplitude spectrum *ρ*(*ω*) and the phase spectrum *θ*(*ω*) obtained from *E*(*t*) with the use of a window function and a zero-filling technique are shown in Fig. 4(b); *ρ*(*ω*) has a broad shape peaking at *ω*/2*π* ~1.5 THz, while *θ*(*ω*) has an almost linear component presumably due to a time-origin misplacement.

The causality-based functions *K*(*δt*) computed with (*ω*_{0}/2*π*, *ω*_{1}/2*π*, *ω*_{2}/2*π*) = (0.47,0.70, 3.16), (0.47, 1.55, 3.16), and (0.47, 2.80, 3.16) THz in an integration frequency range of *ω'*/2*π* = 0.2–5.0 THz are shown in Fig. 4(c) by curves 1–3, respectively. The deviation in *K*(*δt*) from the baseline is found to start at *δt'* = –1.46 ± 0.06 ps [see also the magnified view in Fig. 4(d)]. The actual misplacement is thus estimated to be *δt* = *δt'* – *τ*_{res} = –1.96 ± 0.06 ps, which requires us to correct the time origin to the point indicated by the arrow in Fig. 4(a). Note that the time origin determined here is very close to the position of the largest negative peak and is clearly way from those of the small positive peaks. This result strongly supports an interpretation of the THz emission described in [15]: the THz waveform shown in Fig. 4(a) should have appeared when it was emitted by the creation of a delta-function-like polarization and measured with slight distortion through the electro-optic sensor.

## 4. Summary

We presented a generalized causality-based function with respect to a possible misplacement *δt* of the time origin in THz emission spectroscopy, by deriving a multiply subtractive K-K relation for complex Fourier spectra of emitted THz electric fields. The function was found to have good properties that allow us to determine the time origin in trial THz waveforms simulated for a damped oscillation current and for delta-function-like instantaneous polarizations (created typically in nonlinear optical processes). Furthermore, we demonstrated how the time-origin misplacement can indeed be detected for a recent example of THz waveforms that were emitted from biased conjugated polymers and measured with slight distortion. The function will serve as a useful tool for inferring the time evolution of true current and polarization on the (sub)picosecond scale and allow for a better understanding of ultrafast phenomena in a variety of optoelectronic materials, such as charge separation in photovoltaic devices.

## Funding

JSPS KAKENHI (Grant Numbers JP15K13334 and JP17H03232); Nagaoka University of Technology Presidential Research Grant.

## References

**1. **R. de L. Kronig, “On the theory of dispersion of X-rays,” J. Opt. Soc. Am. **12**(6), 547–557 (1926). [CrossRef]

**2. **H. A. Kramers, “La diffusion de la lumière par les atomes,” in *Atti del Congresso Internazionale dei Fisici, Como* (Zanichelli, 1927), Vol. 2, pp. 545–557.

**3. **V. Lucarini, J. J. Saarinen, K.-E. Peiponen, and E. M. Vartiainen, *Kramers-Kronig Relations in Optical Materials Research* (Springer, 2005).

**4. **V. Lucarini, Y. Ino, K.-E. Peiponen, and M. Kuwata-Gonokami, “Detection and correction of the misplacement error in terahertz spectroscopy by application of singly subtractive Kramers-Kronig relations,” Phys. Rev. B Condens. Matter Mater. Phys. **72**(12), 125107 (2005). [CrossRef]

**5. **T. Unuma, Y. Ino, K.-E. Peiponen, E. M. Vartiainen, M. Kuwata-Gonokami, and K. Hirakawa, “Causality-based method for determining the time origin in terahertz emission spectroscopy,” Opt. Express **19**(13), 12759–12765 (2011). [CrossRef] [PubMed]

**6. **M. Bernier, F. Garet, J.-L. Coutaz, H. Minamide, and A. Sato, “Accurate characterization of resonant samples in the terahertz regime through a technique combining time-domain spectroscopy and Kramers-Kronig analysis,” IEEE Trans. Terahertz Sci. Technol. **6**(3), 442–450 (2016). [CrossRef]

**7. **H. Son, D.-H. Choi, and G.-S. Park, “Improved thickness estimation of liquid water using Kramers-Kronig relations for determination of precise optical parameters in terahertz transmission spectroscopy,” Opt. Express **25**(4), 4509–4518 (2017). [CrossRef] [PubMed]

**8. **T. Unuma, Y. Ino, M. Kuwata-Gonokami, G. Bastard, and K. Hirakawa, “Transient Bloch oscillation with the symmetry-governed phase in semiconductor superlattices,” Phys. Rev. B Condens. Matter Mater. Phys. **81**(12), 125329 (2010). [CrossRef]

**9. **A. Naka, K. Hirakawa, and T. Unuma, “Capacitive response and room-temperature terahertz gain of a Wannier-Stark ladder system in GaAs-based superlattices,” Appl. Phys. Express **9**(11), 112101 (2016). [CrossRef]

**10. **K. Sakai and M. Tani, “Introduction to terahertz pulses,” in *Terahertz Optoelectronics* ed. K. Sakai (Springer, 2005).

**11. **K. F. Palmer, M. Z. Williams, and B. A. Budde, “Multiply subtractive kramers-kronig analysis of optical data,” Appl. Opt. **37**(13), 2660–2673 (1998). [CrossRef] [PubMed]

**12. **V. Lucarini, J. J. Saarinen, and K.-E. Peiponen, “Multiply subtractive Kramers-Krönig relations for arbitrary-order harmonic generation susceptibilities,” Opt. Commun. **218**(4), 409–414 (2003). [CrossRef]

**13. **T. Unuma and A. Matsuda, “Temperature-dependent spectral linewidths of terahertz Bloch oscillations in biased semiconductor superlattices,” Appl. Phys. Lett. **112**(16), 162107 (2018). [CrossRef]

**14. **T. Unuma, Y. Ino, M. Kuwata-Gonokami, E. M. Vartiainen, K.-E. Peiponen, and K. Hirakawa, “Determination of the time origin by the maximum entropy method in time-domain terahertz emission spectroscopy,” Opt. Express **18**(15), 15853–15858 (2010). [CrossRef] [PubMed]

**15. **T. Unuma, N. Yamada, and H. Kishida, “Terahertz emission from biased conjugated polymers excited by femtosecond laser pulses,” Appl. Phys. Express **9**(12), 121601 (2016). [CrossRef]

**16. **Q. Wu and X.-C. Zhang, “7 terahertz broadband GaP electro-optic sensor,” Appl. Phys. Lett. **70**(14), 1784–1786 (1997). [CrossRef]

**17. **G. Gallot, J. Zhang, R. W. McGowan, T.-I. Jeon, and D. Grischkowsky, “Measurements of the THz absorption and dispersion of ZnTe and their relevance to the electro-optic detection of THz radiation,” Appl. Phys. Lett. **74**(23), 3450–3452 (1999). [CrossRef]