Approximations used in the analysis of signals in pump-probe spectroscopy

TL;DR

This paper evaluates the validity of common approximations in pump-probe spectroscopy, especially optical pump-THz probe, revealing they are only accurate in extreme optical thickness cases.

Contribution

It introduces a general transfer matrix method to assess when typical approximations in pump-probe spectroscopy are valid.

Findings

Approximations are valid only in extreme optical thickness regimes.

The transfer matrix approach provides a systematic evaluation.

Most common approximations are not generally applicable.

Abstract

In pump-probe spectroscopy, one often needs to analyse the transmission or reflection of electromagnetic waves through optically pumped media. Here, it is common practice to approximate the analysis in order to extract non-equilibrium values for the optical constants of the media. Concentrating on optical pump-THz probe spectroscopy, and using the transfer matrix approach, we present a general method for evaluating the applicability of the most common approximations. Somewhat surprisingly, we find that these approximations are truly valid only in extreme cases where the optical thickness of the sample is several orders of magnitude smaller or larger than the probe wavelength.

Equations

$$\left(\begin{array}{c}{E}_{j-1,j}\\ Z_0{H}_{j-1,j}\end{array}\right)=\left(\begin{array}{cc}\cos {\delta }_j& \frac{-i\sin {\delta }_j}{{\gamma }_j}\\ -i{\gamma }_j\sin {\delta }_j& \cos {\delta }_j\end{array}\right)\left(\begin{array}{c}{E}_{j,j+1}\\ Z_0{H}_{j,j+1}\end{array}\right),$$

$$M_j\equiv \left(\begin{array}{cc}\cos {\delta }_j& \frac{-i\sin {\delta }_j}{{\gamma }_j}\\ -i{\gamma }_j\sin {\delta }_j& \cos {\delta }_j\end{array}\right).$$

$$\left(\begin{array}{c}{E}_{12}\\ Z_0{H}_{12}\end{array}\right)$$

$$=M_{tot}\left(\begin{array}{c}{E}_{N,N+1}\\ Z_0{H}_{N,N+1}\end{array}\right),$$

$$M_{tot}=\left(\begin{array}{cc}{m}_{11}& m_{12}\\ m_{21}& m_{22}\end{array}\right).$$

$$\left(\begin{array}{c}{E}_1^i+E_1^r\\ (E_1^i-E_1^r){\gamma }_1\end{array}\right)=M_{tot}\left(\begin{array}{c}{E}_{N+1}^t\\ {\gamma }_j{E}_{N+1}^t\end{array}\right).$$

$$r$$

$$t$$

$$\frac{E_{exc}(\omega ,{ϵ}_{exc})}{E(\omega )}=\frac{t_{exc}(\omega ,{ϵ}_{exc})}{t(\omega )},$$

$$\frac{\Delta E(\omega ,{ϵ}_{exc})}{E(\omega )}=\frac{t_{exc}(\omega ,{ϵ}_{exc})}{t(\omega )}-1.$$

$$\Delta ϵ(z)$$

$$\Delta {ϵ}_{short}$$

$$\Delta {ϵ}_{long}$$

$$n^2(\omega )=ϵ(\omega )=1+\frac{i\sigma (\omega )}{\omega {ϵ}_0},$$

$$t_{ij}$$

$$r_{ij}$$

$$E$$

$$E_{exc}$$

$$\frac{\Delta E}{E}$$

$$=\frac{t_{1e}{t}_{e3}{e}^{{\beta }_{\Delta }}}{t_{12}{t}_{23}}\frac{(1-r_{21}{r}_{23}{e}^{2{\beta }_2})}{(1-r_{e1}{r}_{e3}{e}^{2({\beta }_2+{\beta }_{\Delta })})}-1,$$

$$\frac{\Delta E}{E}=$$

$$\times \frac{1-r_{21}{r}_{23}(1+2{\beta }_2)}{1-r_{e1}{r}_{e3}(1+2{\beta }_2+2{\beta }_{\Delta })}-1$$

$$\Delta n=\frac{c(c(n_1+n_3)+id(n_1-n_2)(n_2-n_3)\omega )}{d\omega (2ic{n}_2+d(n_1-n_2)(n_3-n_2)\omega )}\frac{\Delta E}{E}$$

$$\Delta n$$

$$\frac{\Delta E}{E}$$

$$\frac{t_{1e}{t}_{e1}}{t_{12}{t}_{21}}$$

$$\frac{(1-(r_{21}{e}^{{\beta }_2}{)}^2)}{(1-(r_{e1}{e}^{{\beta }_2}{)}^2(1+2{\beta }_{\Delta }))}\approx 1-MR\times \Delta n$$

$$MR=\frac{2i{e}^{2{\beta }_2}(n_2-n_1)(-2ic{n}_1+d\omega (n_2^2-n_1^2))}{c{e}^{2{\beta }_2}(n_1-n_2{)}^2(n_1+n_2)-c(n_1+n_2{)}^3}$$

$$\Delta n={\left[\frac{i\omega d}{c}+\frac{n_2-n_1}{n_2(n_2+n_1)}+MR\right]}^{-1}\frac{\Delta E}{E},$$

$$\Delta n$$

$$\Delta {ϵ}_{short}$$

$$\Delta {ϵ}_{long}$$

$$\Delta {\sigma }_{short}$$

Full text

††journal: josab

Optical pump - THz probe time-domain spectroscopy (OPTP) is a powerful tool for studying the ultrafast charge-carrier dynamics of materials on a femtosecond to nanosecond time-scale, and has successfully been applied to a range of different materials including semiconductors, carbon nanotubes, nanowires, and polymers [1, 2, 3, 4]. The power of this method lies in the extraction of the complex material parameters (represented by a complex permittivity, refractive index or conductivity) from photoinduced changes to the THz transmission/reflection of a sample [5, 6, 7]. This approach has allowed the study of a wealth of photo-species and phenomena, including exciton dynamics [8], plasmon formation [9], coulomb screening [10] and carrier scattering [11].

Extracting the exact, non-equilibrium material parameters of a sample in an OPTP experiment is not trivial, requiring full wave solutions to Maxwell’s equations, which must be numerically solved [12, 13]. However, depending on the geometry and dimensions of the sample, in some cases a set of simplified equations may be used to obtain approximate material parameters of the photoexcited sample, as shown in [7]. The most commonly used approaches assume either a sample which is thin or thick relative to the wavelength of the incident THz radiation, and have been used by numerous groups throughout the literature for a variety of thin- [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 12, 32, 33, 34, 6, 35, 36, 37, 38, 3, 39, 40, 41, 42, 43, 44, 45, 7, 10, 46, 47] and thick-films [48, 49, 6, 50, 51, 52, 53, 1]. Less simplified approximations have also been used for more complex geometries [3, 54, 8, 55, 56, 57, 58]. All of these approximations are based on a number of assumptions and approximations in terms of the probe wavelength and the magnitude of the relative change in transmission/reflection due to photoexcitation. It is important to note that improper use of these simplified equations can give a significant deviation from the expected material parameters, which may in turn lead to incorrect conclusions regarding the underlying charge-carrier dynamics of the material [6, 59, 12, 60]. However, assessing the appropriate approximation for a given sample is not always trivial, since there are no clear boundaries for validity. Despite this uncertainty, only a few groups have investigated the validity of these approximations in depth, and only for few specific cases in the reflection geometry [59, 12]. Using a 0.33 mm undoped GaAs wafer as a case study, D’Angelo et al. [59] demonstrated that the approximations of a step-like excitation and a thin film approximation are generally insufficient for extracting the transient optical properties of thin films in the reflection geometry, due to the strong phase-sensitivity of these measurements. Similarly, Hempel et al. [12] investigated the effect of substrate refractive index and thickness on the sensitivity of thin film transient reflection measurements. By comparing thin films of polycrystalline Cu2SnZnSe4 placed on a glass substrate and on a highly conductive substrate of molybdenum, they showed that the thin film approximation breaks down when film and substrate refractive index differs significantly. While these studies highlight potential problems with specific approximations and sample geometries, it can be difficult to determine to what extent problems persist in other, more common sample geometries.

In this paper, we present a general method for evaluating the accuracy of various approximations for any particular sample geometry. It should be noted that regimes of applicability can be evaluated using series expansions of the transmission functions for the given geometry, since many of the commonly used approximation rely on first order Taylor expansions. However, this is not straight forward even for the simple case of a single homogeneously photoexcited layer, see supplementary material section S2. Instead, using a modified transfer matrix method, we simulate the transmitted and reflected fields of a photoexcited, multilayer structure, where the thickness of each layer and the photoinduced change in the complex permittivity are treated as input parameters. We then use these fields to extract approximate solutions for the complex permittivity of the material, which we compare to the input parameters, allowing us to assess ranges of validity. Using this method we calculate the relative deviation of the most commonly used approximations for a broad range of sample geometries, where the probe wavelength ranges from much smaller to much larger than the optical penetration depth of the pump beam in the material. Somewhat surprisingly, we find that these approximations are truly valid only in extreme cases where the optical thickness of the sample is several orders of magnitude smaller or larger than the probe wavelength. We then go beyond the most commonly applied approximations, using a general numerical approach based on the transfer matrix method, which exhibits a much broader range of validity.

1 The Transfer Matrix Method

In a standard OPTP experiment the material of interest will typically be part of some multilayer structure. For example, the sample can be placed on a substrate or in a cuvette. On photoexcitation, part or all of this multilayer structure experiences a change in optical parameters, and the photoexcited change in the THz transmission/reflection of the entire structure is measured. Note that in this article we deal only with cases where the photoexcited changes to the optical parameters vary on a timescale which is slow compared to the period of the THz field, where the excited state optical parameters can be treated as quasi-static. In such circumstances, one wants a full wave solution capable of determining the complex transmission and reflection coefficients of arbitrary multilayers of differing index. The most straightforward approach here uses the transfer matrix method (TMM). Below, we summarize the important definitions and equations, but a full derivation can be found in many introductory textbooks on optics, for example Pedrotti et al. [61, Chapter 22,].

Figure 1 shows a general $N$-layer structure consisting of homogeneous layers of thickness $d_{j}$ and complex permittivity $\epsilon_{j}$, where $j=1,2,3,\dots,N+1$. Here layer 1 and $N+1$ are the semi-infinite incident and transmitted regions, respectively. The electric ($E$) and magnetic ($H$) fields at interfaces $j$ and $j+1$ are related through

[TABLE]

where $Z_{0}$ is the free-space impedance, $\delta_{j}=\omega\sqrt{\epsilon_{j}}d_{j}\cos{\theta_{j}}/c$ is the complex phase (including absorption) accumulated through the $j$’th layer, $\theta_{j}$ is the incident angle on the $j$’th interface, $\omega$ is the angular frequency of the fields, $c=299792458$ m/s is the vacuum speed of light, and $\gamma^{TE}_{j}=\sqrt{\epsilon_{j}}\cos{\theta_{j}}$ and $\gamma^{TM}_{j}=\cos{\theta_{j}}/\sqrt{\epsilon_{j}}$ for TE- and TM-polarization, respectively. We define the transfer matrix for the $j$-th layer as

[TABLE]

Note that, in the absence of gain, since the off axis components of the transfer matrix have a negative sign, the formulation above is consistent with a positive sign for the imaginary part of the permittivity as per definition. The fields transmitted through the entire multilayer structure are related to the incident fields via the relation

[TABLE]

where $M_{tot}$ is the total transfer matrix of the system, given by

[TABLE]

(3) can be rewritten in terms of $E^{i}_{1}$, $E^{r}_{1}$ and $E^{t}_{N+1}$:

[TABLE]

The reflection ($r$) and transmission ($t$) coefficients for the entire structure can then be found from (5)

[TABLE]

If $\epsilon_{j}$ and $d_{j}$ is known for each layer, it is fairly straight-forward to calculate the transmitted and reflected fields of the multilayer structure. The reverse is also possible, i.e. (6) or (7) can be solved for an unknown $\epsilon_{j}$, assuming the incident and transmitted/reflected fields are known, which is normally the case in OPTP experiments. However, without approximations to (6) and (7) their solutions must be found numerically, which can be resource and time intensive, depending on the complexity of the multilayer structure.

2 Wave-propagation in an optically pumped medium

The goal of an OPTP experiment is normally to determine the photoinduced change in the THz optical properties of a sample. This is done by measuring either the transmitted or reflected electric field of the photoexcited sample $E_{exc}$, as well as an appropriate reference, usually the unexcited sample $E$. From these two measurements, the non-equilibrium permittivity of the photoexcited sample $\epsilon_{exc}$ can be determined by solving the equation

[TABLE]

where $t$ and $t_{exc}$ are the complex transmission functions of the unexcited and photoexcited multilayer structure. (8) is commonly written in terms of the change in transmitted electric field $\Delta{}E=E_{exc}-E$:

[TABLE]

Assuming the equilibrium permittivity and thickness of each layer of the multilayer structure are known, and that the dynamics of relaxation are slow compared the period of the THz radiation, one can solve (9) as a function of the non-equilibrium permittivity $\epsilon_{exc}$ for a given angular frequency $\omega$. In figure 2 we consider a common scenario of a sample with permittivity $\epsilon_{2}$, surrounded by an incidence ($\epsilon_{1}$) and transmission medium ($\epsilon_{3}$). If the sample is photoexcited, one expects a non-equilibrium permittivity that varies spatially through the sample, depending on the penetration of the pump pulse into the sample. If pump absorption is linear, and the amplitude of the THz permittivity also varies linearly with the density of photospecies, the photoinduced change in permittivity $\Delta\epsilon$ is expected to decay exponentially in the propagation direction, following the attenuation of the incident pump, and described by a decay length defined by a penetration depth $d_{p}$ of the pump light, as assumed in [7, 62, 1]:

[TABLE]

where $\Delta\epsilon_{s}$ is the change in permittivity at the surface of the sample, which is assumed to be sufficiently small so that it scales linearly with the incident pump light. In this case, the spatially varying permittivity can be approximated by dividing the sample into $N$ homogeneous layers, each with a permittivity $\epsilon_{exc}(z)$ determined by the distance $z$ into the sample. Using (7) and (9) and (10) it is then possible to determine $\Delta{}\epsilon_{s}$ from the experimentally obtained $\Delta{}E/E$, given that thickness ($d$), penetration depth ($d_{p}$) and equilibrium permittivity ($\epsilon_{2}$) of the sample is known. However, no analytical solution exists for $\Delta{}\epsilon_{s}$ in its current form, so a solution must be obtained numerically, which can be resource and time consuming depending on the size of $N$.

A common approximation is to represent the photoexcitation as a single homogeneous layer of thickness $d_{p}$ and constant photoexcited permittivity $\epsilon_{exc}$, i.e. a step-like excitation approximation. In this case, depending on the size of the wavelength inside the sample $\lambda_{s}=\lambda_{0}/Re(\sqrt{\epsilon_{2}})$ compared to the thickness and penetration depth, various further approximations can be made, which result in a simpler and analytical solution for $\Delta{}\epsilon$ [7, 6, 51]. In order to obtain an analytical solution, the most common approximations assume a small perturbation of the photoexcited sample, i.e. $\Delta{}E/E<<1$, and that the wavelength is much shorter or much larger than the thickness and penetration depth of the sample, i.e. $|\beta|^{-1}=\frac{\lambda_{0}}{2\pi d|\sqrt{\epsilon_{2}}|}<<1$ [51, 52, 53] or $\frac{\lambda_{0}}{2\pi d|\sqrt{\epsilon_{2}}|}>>1$ [6, 35, 36, 26, 37, 38, 3, 39, 40, 41, 42, 28, 43, 44, 45], resulting in

[TABLE]

(11) is known as the short wavelength limit or thick sample approximation, and (12) is known as the long wavelength limit or thin-film approximation. A full derivation can be found in appendix S1, along with equivalent expressions written in terms of either the complex refractive index or the complex conductivity. If Re$(\sqrt{\epsilon_{2}})>$ Im$(\sqrt{\epsilon_{2}})$, then the conditions for the short and long wavelength limits approximately become $\lambda_{s}/d<<1$ and $\lambda_{s}/d>>1$, respectively. However, as we will show later on, $\lambda_{s}/d$ is actually not a good parameter for judging the validity of a given approximation and can be a bit misleading. Note that (11) and (12) assume a homogeneous excitation of the entire sample, which occurs for weakly absorbing materials, i.e. $d_{p}\approx d$. For strongly absorbing materials, i.e. $d_{p}<<d$, similar expressions to (11) [48, 49, 1] and (12) [14, 25, 18, 20, 29, 7, 10, 27, 46, 47] can be derived, where $d$ is replaced by $d_{p}$. In this case, the photoexcited region is approximated as thin layer of thickness $d_{p}$, located at the surface of the sample.

We now examine the relative deviation of the four different analysis approaches mentioned so far. Specifically: (i) the short wavelength limit of (11) with $d=d_{p}$; (ii) the long wavelength limit of (12) with $d=d_{p}$; (iii) the TMM approach with the photoexcited region of the sample approximated as a single homogeneous layer of thickness $d_{p}$ and permittivity $\epsilon_{exc}$; and (iv) the full TMM approach using $N=10$ layers to represent the photoexcited sample. In order to analyse the relative deviation of the different approximations, we compare results for a sample of known thickness $d$, penetration depth $d_{p}=d/4$, and equilibrium permittivity $\epsilon_{2}=10.7+6.8i$, with the incident region being air ($\epsilon_{1}=1$), and the transmitted region being a quartz substrate ($\epsilon_{3}=3.8$) [63]. We assume a photoexcited sample response described by $\Delta{}\epsilon(z)$, and modelled according to (10), with $\Delta{}\epsilon_{s}=(1+1i)\times 10^{-6}$. Note that such a sample is representative of a semiconductor thin-film deposited on a THz-transparent substrate. This type of sample has often been examined using terahertz spectroscopy, for example in the case of semiconductor nanowires [3]. As a simple example we have chosen the value of $\epsilon_{2}$ to be that of doped silicon [64] at 1 THz, since this semiconductor material has a noticeable absorption in the THz frequency range, and we have chosen $\Delta\epsilon_{s}$ to be sufficiently small so that $\Delta{}E/E<<1$ for all wavelengths ranging from $\lambda_{s}/d=10^{-3}$ to $\lambda_{s}/d=10^{5}$. We find that the transmission of the sample is essentially constant for $N>100$, hence we use calculations with $N=100$ layers to be representative of the true transmission. We then feed the generated $\Delta{}E/E$ into the various approximate analyses and compare the resulting $\Delta{}\epsilon_{s}^{out}$ with the input value of $\Delta{}\epsilon_{s}^{in}$. This gives a relative deviation $(\Delta\epsilon_{s}^{out}-\Delta\epsilon_{s}^{in})/\Delta\epsilon_{s}^{in}$ with which to assess the validity of the various approximations described above. Note that the code used for this analysis can be found on the ArXiv supplementary material repository for this article, and can be adapted for other sample geometries and tested using other approximations.

The resulting magnitude of the relative deviation is plotted in figure 3a as a function of the wavelength normalized to the thickness of the sample $\lambda_{s}/d$. Note that in this context, the magnitude of the relative deviation tells us how close $\Delta\epsilon_{s}^{out}$ is to $\Delta\epsilon_{s}^{in}$, while the phase tells us whether $\Delta\epsilon_{s}^{out}$ is smaller or larger than $\Delta\epsilon_{s}^{in}$ and if the ratio of the real and imaginary parts are distributed correctly. A corresponding phase plot of the data in figure 3 can be seen in figure S2 in the supplementary material. As one might expect, there is a significant error when the wavelength is comparable with the sample thickness, i.e. $\lambda_{s}/d\approx 1$, with up to 20–50% in the magnitude of the relative deviation for the approximations in (11) and (12), and 11% relative deviation when using the TMM approach assuming a homogeneous photoexcited layer. The short and long wavelength approximations converge to a 2% relative deviation for $\lambda_{s}/d<<1$ and $\lambda_{s}/d>>1$, respectively. This asymptotic value corresponds more or less exactly with the factor $1-(1-exp[-d/d_{p}])=exp[-d/d_{p}]=1.8\%$ found in eq. (27) from ref. [7], which the simple long- and short-wavelength limits do not account for. Indeed, when varying $d/d_{p}$, we observe that the asymptotic values shift accordingly, see supplementary material section S3 for more details. However, more surprisingly, the relative deviation of both the short and long wavelength approximations is considerably larger than we expected in the region $\lambda_{s}/d$ = $10^{-1}$ to $10^{3}$. For example, the relative deviation of the long wavelength limit, even when the wavelength is two orders larger than the sample thickness ($\lambda_{s}/d$ = $10^{2}$), is 10 %. This underlines the problem in using these approximations without properly analysing their applicability for a given sample geometry, as has previously been observed by D’Angelo et al. [59] and Hempel et al. [12]. Furthermore, we note that the error of these approximations will scale with the relative variation between $\epsilon_{2}$ and the surrounding mediums ($\epsilon_{1},\epsilon_{3}$), in agreement with similar observations by Hempel et al. [12], who showed that a high variation between sample and substrate refractive indices can lead to significant errors in the long wavelength approximation in the reflection geometry. To avoid this issue, one can apply the TMM approach in full to better approximate the spatially varying profile of the non-equilibrium optical constants, as can be in figure 3a, where for $N=10$ we obtain a relative deviation of less than 2% at $\lambda_{s}/d=1$. In figures 3b-3c we plot a similar analysis carried out using the TMM approach while varying the number of layers used to approximate the exponentially varying permittivity. We see that for $\lambda_{s}/d=1$, the TMM approach converges relatively quickly for $N<5$ in terms of the magnitude of the relative deviation, but slows down beyond that, achieving a relative deviation less than of 1% for $N\geq 12$. This demonstrates that the TMM approach is valid for all wavelengths as long as $N$ is large enough and thus TMM is a valid option in the case where no other approximation is applicable. However, we note that there is a trade off in terms of the computation time and memory requirement for this approach, which will scale non-linearly with $N$. Furthermore, when solving equations numerically there is always a chance that the obtained solution is not the correct one, but rather a local minimum, especially if the initial guess used in these routines is not close to the correct solution. To avoid this we solve (9) numerically in terms of the photoexcited permittivity $\epsilon_{exc}$, using the unexcited permittivity $\epsilon_{2}$ as the initial guess, which is a fair approximation when $\Delta{}E/E<<1$ (and therefore $\Delta\epsilon_{s}/\epsilon_{2}<<1$). However in the far majority of cases, one of the approximations found in ref. [7] can be used to extract the proper photoexcited permittivity of the sample, see supplementary material section S3 where we compare the performance of Eq. (27) from ref. [7] with the long wavelength approximation. Unlike the TMM approach, the approximations found in ref. [7] also account for phase matching and frequency mixing effects, which must be considered when the decay-time of the photoexcitation is similar to the duration of the THz pulse. However in this context, the method presented in this paper can be used as a quick comparison to work out which approximation is the most appropriate for a given sample geometry.

Finally in figure 4 we examine a specific case of experimentally obtained data from a typical OPTP experiment [1], in this case a film of single-walled carbon nanotubes of thickness 700 nm on a quartz substrate. We take $\Delta{}E/E$ measured two picoseconds after photoexcitation with an 800 nm, 100 fs pump pulse and incident fluence of 52 $\mu$J/cm2, see [65] for details of the sample preparation and the experimental setup. The sample geometry is the same as in figure 2, where the permittivity of the sample at 1 THz is approximately $\epsilon_{2}=88+234i$, $\epsilon_{1}=1$ and $\epsilon_{3}=3.8$. In this case $\Delta{}E/E\leq 0.5\%$, $\lambda_{s}/d=37$ for 1 THz, and the penetration depth $d_{p}$ is 285 nm. Observing the sample geometry (a thin film of thickness ten times thinner than the wavelength), one might assume the long wavelength limit in (12) to be a valid approximation in this case. In figure 4 we plot the obtained complex photo-induced change in permittivity $\Delta{}\epsilon_{s}$, extracted using the long wavelength limit, which we compare with the correct non-equilibrium permittivity $\Delta{}\epsilon_{s}$ obtained using the full TMM approach with $N=10$, which we ensured had converged for this valued of $N$, see supplementary material section S4 for details. We see that the two methods give remarkably different spectral signals and magnitudes, illustrating the problems that can occur with improper use of approximate analysis methods. In this case, neither of the conditions $d_{p}/d<<1$ or $\lambda_{0}/(2\pi d\sqrt{\epsilon_{2}})>>1$ are properly satisfied for use of (12). Such a difference in the sample response could well result in misinterpretation of the photo-response, as discussed in ref. [6]. In this particular case, the change in sign of $\Delta{}\epsilon_{s}$ observed at 0.4 THz in figure 4b is indicative of a plasmon resonance in this material [66, 65, 67], the parameters of which would be badly misinterpreted if one applied an approximate analysis.

3 Conclusion

In conclusion we have formulated a general method for evaluating the applicability of the most common approximation used in optical pump-THz probe spectroscopy. Somewhat surprisingly, we find that these approximations are truly valid only in extreme cases where the optical thickness of the sample is several orders of magnitude smaller or larger than the probe wavelength.This demonstrates the need to properly analyse the applicability of an approximation for a given sample geometry, since improper use of these approximations can lead to significant errors, which in turn can lead to an incorrect interpretation of the photoexcited response.

Acknowledgements

This research was partially supported by the European Union’s Seventh Framework Programme (FP7) for research, technological development and demonstration under project 607521 NOTEDEV. The Mathematica and Matlab code used for performing the Transfer Matrix calculations was initially adapted from the Mathematica code written by Steven Byrnes [68]. The Mathematica code used for this analysis can be found on the ArXiv supplementary material repository for this article.

Appendix S1 Derivation of the long- and short-wavelength limit

The short and long wavelength limits in (11) and (12) can be derived from the Fresnel equations under the assumption that the photo-induced change in the dielectric properties of the sample are small, i.e. $\Delta{}E/E<<1$. Furthermore, we assume a homogeneous excitation of the entire sample shown in figure 2. When discussing the dielectric properties of a sample, it is common to consider either the complex refractive index, $n$, or the complex conductivity, $\sigma$, which are related to the permittivity, $\epsilon$, through

[TABLE]

where $\epsilon_{0}=8.85\times 10^{-12}Fm^{-1}$ is the vacuum permittivity and $\omega$ is the angular frequency. For the purpose of deriving (11) and (12), we use the refractive index $n$, since it simplifies the equations.

The Fresnel transmission and reflection coefficients for normal incidence are given by

[TABLE]

For the geometry shown in figure 2 the transmitted electric field through the photoexcited ($E_{exc}$) and unexcited ($E$) sample is then given by

[TABLE]

where the indices $e$ and $s$ denotes the photoexcited and unexcited sample, respectively, and $\beta_{j}=i\omega n_{j}d/c=i\delta_{j}$ is the accumulated phase through layer $j$. In (16) and (17), the numerator represents direct transmission through the sample, while the denominator represents multiple internal reflections in the sample. The relative change in the transmitted electric field $\Delta{}E/E$ is then given by

[TABLE]

where $\beta_{\Delta}=i\omega\Delta{n}d/c$ and $\Delta{}n=n_{exc}-n_{2}$.

(18) cannot be solved analytically for $\Delta{}n$ in it’s current form, however we can approximate the equation in several ways, depending on whether $\omega|n_{2}|d/c<<1$ (long wavelength limit) or $\omega|n_{2}|d/c>>1$ (short wavelength limit).

In case of the long wavelength limit, we can Taylor-expand all the exponentials in (19), since $|\beta_{2}|,|\beta_{\Delta}|<<1$:

[TABLE]

(S1) can then be reduced by discarding higher order terms of $\Delta{}n$, and subsequently solved for $\Delta{}n$:

[TABLE]

By substituting $x=\omega d/c$ in (21) and Taylor-expanding in $x$ up to the first order, the equivalent approximation of (12) is obtained in terms of $\Delta{}n$:

[TABLE]

for $\omega d\rightarrow 0$.

In case of the short wavelength limit, the assumption $|\beta_{2}|<<1$ no longer holds, so a slighty different approach is required, which is outlined in ref. [51] and repeated in this appendix. In order to simplify the derivation, we begin by assuming $n_{3}=n_{1}$, i.e. the incident and transmitted medium are the same. It is possible to derive (11) without this assumption, however the derivation becomes much more cumbersome while the final equation will be the same in either case. Using $n_{3}=n_{1}$, (19) becomes

[TABLE]

To approximate (23), we must rely on the assumption $\Delta{}n/n<<1$. We begin by Taylor-expanding the exponential $e^{\beta_{\Delta}}$, since the assumption $\beta_{\Delta}<<1$ is still valid:

[TABLE]

Similarly, the transmission and the multiple reflection terms can be approximated by Taylor expanding in terms of $\Delta{}n$:

[TABLE]

and

[TABLE]

where the multiple reflection term $MR$ is given by

[TABLE]

From equations (24), (25) and (26) we can then obtain an analytical solution for $\Delta{}n$:

[TABLE]

where the first term represents modification of the propagation through the sample, the second term represents changes in reflective losses at the two interfaces, and the third term represents the contribution from multiple reflections. For very thick samples, i.e. the short wavelength limit ($\omega d/c>>1$), the propagation term dominates (28) and thus the other terms can be disgarded, giving us the equivalent short wavelength limit of (11) written in terms of $n$:

[TABLE]

for $\omega d\rightarrow\infty$.

(29) and (22) can be rewritten in terms of the complex permittivity using $\Delta\epsilon=(n_{2}+\Delta{}n)^{2}-n_{2}^{2}$, resulting in (11) and (12):

[TABLE]

where we have discarded the second order term of $\Delta{}E/E$ under the assumption that $\Delta{}E/E<<1$. Likewise, the complex conductivity can be found from (30) and (31) using $\Delta\sigma=-i\omega\epsilon_{0}\Delta\epsilon$:

[TABLE]

It is important to note that for the long wavelength limit in (22), (31) and (33) we have assumed that the wavelength is much larger than the thickness of the sample. However, it is straightforward to derive a similar expression, where we consider a thin photoexcited region of thickness $d<<\lambda$ at the surface of a semi-infinite sample, which is another common approximation in OPTP. In this case, $d$ becomes the thickness of the photoexcited region $d_{p}$. Likewise the $d$ in (29), (30) and (32) can be interchanged with the penetration depth $d_{p}$ in the case where we are no longer considering a homogeneous excitation of the entire sample.

Appendix S2 Assessment of Approximations through Taylor Expansion

In this section we demonstrate an alternative method for assessing a given approximation by considering the higher order terms of the Taylor expansion of the transmission function for a given sample geometry. As an example we consider the simplest case of a homogeneously photoexcited sample in air, which is described by (19):

[TABLE]

We established previously that for $|\beta_{2}|,|\beta_{\Delta}|<<1$, the long wavelength limit (12) is a valid approximation. Since this approximation relies on Taylor expanding the exponentials in (34) to the first order, we can assess the validity of this approximation by Taylor expanding (34) in terms of $x=\omega d/c$:

[TABLE]

where we have refrained from writing the higher order terms explicitly here due to the complexity of these terms. Using the same input parameters as in the main manuscript (except $d_{p}>>d$ in this case), we then compare the different order terms in (35) with the relative deviation of the long wavelength limit in figure S1:

We observe that the point at which the higher order terms start to dominate in figure S1b is more or less the same point at which the relative deviation of the long wavelength limit becomes massive. The benefit of this approach is that the dependence on the various sample parameters become much clearer, compared to the method presented in the main manuscript, however this method requires clear knowledge of the assumptions and derivation for each approximation, as well as writing out the Taylor expansion explicitly for transmission function of the sample geometry of interest, which can very quickly become tedious and massive for more complex geometries such as the one presented in figure 2.

Appendix S3 Additional Approximations and Tests

In figure S2a we plot the phase of the relative deviation for the same data presented in figure 3, along with the convergence tests in figures S2b-S2c.

As mentioned in the main manuscript, it should be noted that in this context, the magnitude of the relative deviation tells us how close $\Delta\epsilon_{s}^{out}$ is to $\Delta\epsilon_{s}^{in}$, while the phase tells us whether $\Delta\epsilon_{s}^{out}$ is smaller or larger than $\Delta\epsilon_{s}^{in}$ and if the ratio of the real and imaginary parts are distributed correctly. For example, if $\Delta\epsilon_{s}^{out}<\Delta\epsilon_{s}^{in}$, but the ratio of their real and imaginary parts are the same, then the phase of their relative deviation will be $\pm\pi$.

To verify the asymptotic behaviour of the short and long wavelength limits in figure 3, we compare the same approximations for different values of $d/d_{p}$ in figure S3:

We observe that the asymptotic value of the two approximations fit more or less exactly with the factor exp($-d/d_{p}$) for the different values of $d/d_{p}$, which is present in Eq. (27) and similar equations from ref. [7] that accounts for when $d\approx d_{p}$. This indicates that the approximations presented in [7] are generally more accurate widely applicable, which we demonstrate in figure S4:

Here Eq. (27) from ref. [7] can be written in terms of $\Delta\epsilon$ as

[TABLE]

where we have ignored frequency-mixing and phase-matching effects. Interestingly, while (36) is more accurate for large $\lambda_{s}/d>10^{2}$, below this value the two approximations give the same result, indicating that there is still a large wavelength range for which the more complex (36) also becomes unreliable.

Appendix S4 Convergence Tests

To ensure the validity of our TMM simulation, we tested how our results converge as we increase the number of layers $N_{sim}$ used when initially simulating the transmitted fields, see figure S5. We observe that the results fairly quickly converge for $N_{sim}$ greater than approximately 40-60 layers, and become more or less indistinguishable when $N_{sim}$ approaches 100 layers.

To check the validity of our experimental data, we performed a similar convergence test with respect to the number of layers $N$ used in our TMM approach to extract our experimental data in figure 4:

As can be seen in figure S6, the data points for $N\geq 5$ are practically indistinguishable, which fits with our expected relative deviation of a few percent seen in figure 3b. This means that the TMM solution shown in figure 4 is more or less "exact" in this context.