Because the math may be somewhat dense, I outline the procedure below:
- Assume we have an input field with a known spatio-temporal distribution.
- Fourier transform that input field in space and time to get the
plane wave components.
- Apply Fresnel analysis to each component and obtain the plane wave
spectrum for the evanescent side of the interface.
- Superpose all the evanescent components back together (if possible).
- The superposition should result in a single penetration depth, which
may in fact be infinite if all components of the input field are not
above the critical angle.
- Find the resulting exponential decay term and take $\frac{1}{e}$
of that value to be the penetration depth.
Suppose that we have some incident electric field incident on an interface completely polarized in the p-plane, given by $f({\bf{r}},t){\hat{\bf{p}}}$, where ${\hat{\bf{p}}}$ is the unit vector in the direction perpendicular to the Poynting vector for the propagating fields. We also assume that the incident field contains no evanescent wave components. Any vector field whose components are elements of the Schwartz class (of functions) can be represented as a superposition of plane waves \cite{Lax02}, so $f({\bf{r}},t)$ can represented as
\begin{align}
f({\bf{r}},t)=\int_{-\infty}^\infty F({\bf{k}},\omega)e^{i({\bf{k}}\cdot{\bf{r}}-\omega t)} d^3kd\omega
\end{align}
where $F({\bf{k}},\omega)$ is the 4-dimensional Fourier transform of $f({\bf{r}},t)$, (3 spatial dimensions and 1 temporal dimension). In general, ${\bf{k}}$ can be a complex valued vector, but here we assume that the ${\bf{k}}$ is propagating, (i.e., homogeneous) so it is real valued. Additionally, because of the symmetries of a planar interface, we can assume $k_y=0$, and $k_z\geq 0$, when we define the interface to lie in the $xy$-plane, with the Poynting vector of the incoming beam point somewhere in the $+z$ direction, and when, for simplicity, we assume that the beam is only varying spatially in the $xz$-plane. For this case, the general Fresnel amplitude coefficient for the p-polarization is given by :
\begin{align}
\tau_{\text{p}}(k_x,\omega) = &2\frac{\epsilon_a(\omega)\sqrt{\epsilon_b(\omega)\mu_b(\omega) - \left(\frac{ck_x}{\omega}\right)^2}}{\epsilon_a(\omega)\sqrt{\epsilon_b(\omega)\mu_b(\omega)-\left(\frac{ck_x}{\omega}\right)^2}+\epsilon_b(\omega)\sqrt{\epsilon_a(\omega)\mu_a(\omega)-\left(\frac{ck_x}{\omega}\right)^2}}\\
=&2\frac{\epsilon_a(\omega)\sqrt{n_b^2(\omega) - \left(\frac{ck_x}{\omega}\right)^2}}{\epsilon_a(\omega)\sqrt{n_b^2(\omega)-\left(\frac{ck_x}{\omega}\right)^2}+\epsilon_b(\omega)\sqrt{n_a^2(\omega)-\left(\frac{ck_x}{\omega}\right)^2}}\\
\end{align}
where $n(\omega)=\sqrt{\epsilon(\omega)\mu(\omega)}$ is the complex index of refraction, $\epsilon(\omega)$ is the dielectric permittivity, $\mu(\omega)$ is the magnetic permeability, and $\omega$ is angular frequency. Is this case we assume that we are at optical frequencies, i.e., that $\mu(\omega) \approx 1$, then
\begin{align}
n^2(\omega) \to \epsilon(\omega)
\end{align}
which implies that the Fresnel coefficient becomes:
\begin{align}
\tau_{\text{p}}(k_x,\omega) = 2\frac{n_a^2(\omega)\sqrt{n_b^2(\omega) - \left(\frac{ck_x}{\omega}\right)^2}}{n_a^2(\omega)\sqrt{n_b^2(\omega)-\left(\frac{ck_x}{\omega}\right)^2}+n_b^2(\omega)\sqrt{n_a^2(\omega)-\left(\frac{ck_x}{\omega}\right)^2}}\\
\end{align}
From Maxwell's equations we have the following relationship for each plane wave component:
\begin{align}
{\bf{k}}\cdot{\bf{k}}= &\left(\frac{\omega}{c}\right)^2\epsilon(\omega)\mu(\omega)\\
= &\left(\frac{\omega\cdot n(\omega)}{c}\right)^2
\end{align}
where $c$ is the speed of light. In our case this results in:
\begin{align}
k_x^2 + k_z^2= &\left(\frac{\omega \cdot n(\omega)}{c}\right)^2
\end{align}
We now need to answer the question, what is the resultant ${\bf{k}}'$ on the other side of the interface? From the generalized Snell's law we know that, for our case, $k_x = k'_x$. So we must compute $k'_z$ using the generalized Snell's law:
\begin{align}
n_a(\omega)\sin(\theta_a)=&n_b(\omega)\sin(\theta_b)\\
\implies \frac{n_a(\omega)}{n_b(\omega)}\sin(\theta_a)=&\sin(\theta_b)
\end{align}
where $\theta$ is the angle of the ${\bf{k}}$s with respect to the $z$-axis, and the sine functions can be complex valued. If the left hand side of the above is real and $n_a(\omega) > n_b(\omega)$ the we have the following equation where $C$ is real:
\begin{align}
&\frac{e^{i\theta_b}-e^{-i\theta_b}}{2i}=C\\
\implies &e^{i\theta_b}-e^{-i\theta_b}=2Ci\\
\end{align}
In this case we have that $\Re(\theta_b) = \pi /2$, and the imaginary part varies, to meet the above equation. Note that in whatever latex interpreter StackExchange is using $\Re(\cdot)$ is the real part and $\Im(\cdot)$ is the imaginary part. Then
\begin{align}
&ie^{-\Im(\theta_b)}+ie^{\Im(\theta_b)}=2Ci\\
\implies &e^{-\Im(\theta_b)}+e^{\Im(\theta_b)}=2C\\
\implies &\cosh(\Im(\theta_b))=\frac{n_a(\omega)}{n_b(\omega)}\sin(\theta_a)\\
\implies &\Im(\theta_b)=\cosh^{-1}\left(\frac{n_a(\omega)}{n_b(\omega)}\sin(\theta_a)\right)
\end{align}
since $e^{i\pi/2} = i, e^{-i\pi/2}=-1$. Now since $\theta_a$ is defined to the $z$-axis, we have the following definition for ${\bf{k}}$:
\begin{align}
{\bf{k}}=\frac{\omega \cdot n_a(\omega)}{c}\begin{pmatrix}
\sin(\theta_a)\\
0\\
\cos(\theta_a)
\end{pmatrix}.
\end{align}
Note, that we have:
\begin{align}
{\bf{k}}'\cdot {\bf{k}}' = \|{\bf{k}}'_{\text{re}}\| - \|{\bf{k}}'_{\text{im}}\| + 2i{\bf{k}}'_{\text{re}}\cdot {\bf{k}}'_{\text{im}} = \left(\frac{\omega}{c}\right)^2n^2_b(\omega)
\end{align}
for the complex valued ${\bf{k}}'= {\bf{k}}'_{\text{re}} + {\bf{k}}'_{\text{im}}$. Since we are assuming that $n_b(\omega)$ is real valued, then the imaginary part of the dot product must be zero, and hence the real and imaginary parts of ${\bf{k}}'$ are perpendicular. Therefore, since $k'_x = k_x$ and $k_x$ is real, $k'_z$ is either real or purely imaginary. When $k'_z$ is purely imaginary then we have the definition of an evanescent wave. Those components in the input beam which have $k'_z$ purely real will not be evanescent. Now we just need to compute $k'_z$. We know that
\begin{align}
{\bf{k}}'=\frac{\omega \cdot n_b(\omega)}{c}\begin{pmatrix}
\sin(\theta_b)\\
0\\
\cos(\theta_b)
\end{pmatrix}
=\frac{\omega \cdot n_b(\omega)}{c}\begin{pmatrix}
k_x\\
0\\
\cos(\theta_b)
\end{pmatrix}.
\end{align}
It is easy to show that for $\theta_b = \pi/2 + \Im(\theta_b)$ and the above solution for $ \Im(\theta_b)$ that
\begin{align}
\cos(\theta_b) = &i\sinh[\Im(\theta_b)]\\
= &i\sinh\left[\cosh^{-1}\left(\frac{n_a(\omega)}{n_b(\omega)}\sin(\theta_a)\right)\right]\\
= &i\sinh\left[\cosh^{-1}\left(\frac{c}{\omega n_b(\omega)}k_x\right)\right]
\end{align}
Now, we have each ${\bf{k}}$ and associated ${\bf{k}}'$, and the $\tau_p$. The field on the $n_b$ side of the interface is then
\begin{align}
\tau_{\text{p}}(k_x,\omega)F({\bf{k}},\omega)e^{i({\bf{k}}'\cdot{\bf{r}}-\omega t)}
\end{align}
for each ${\bf{k}}$ and each $\omega$. We must superpose these to obtain the field. First, we shall rewrite the above in terms of $k_x$.
\begin{align}
{\bf{k}}'=\frac{\omega \cdot n_b(\omega)}{c}\begin{pmatrix}
k_x\\
0\\
i\sinh\left[\cosh^{-1}\left(\frac{c}{\omega n_b(\omega)}k_x\right)\right]
\end{pmatrix}
\end{align}
Note that the above is only for the evanescent components! We have the following relations (a restatement of above):
\begin{align}
k_x'&=k_x\\
k_y'&=k_y=0\\
k_z'&=i\sinh\left[\cosh^{-1}\left(\frac{c}{\omega n_b(\omega)}k_x\right)\right]
\end{align}
The superposition then becomes:
\begin{align}
\int_{-\infty}^\infty\int_{\bf{E}}\tau_{\text{p}}(k_x,\omega)F(k_x,k_z,\omega)e^{-\sinh\left[\cosh^{-1}\left(\frac{c}{\omega n_b(\omega)}k_x\right)\right]z}e^{i(k_x x-\omega t)}dk_x dk_z dt + \big[\text{propagating components}\big]
\end{align}
where $\bf{E}$ is the domain for the evanescent components. Then by analyzing the inner integral (wrt $k_x$) we get
\begin{align}
&\int_{-\infty}^\infty\mathrm{rect}\left(\frac{(k_x-k_{x,0})}{D}\right)\tau_{\text{p}}(k_x,\omega)F(k_x,k_z,\omega)e^{-\sinh\left[\cosh^{-1}\left(\frac{c}{\omega n_b(\omega)}k_x\right)\right]z}e^{ik_x x}dk_x\\
=&\mbox{$\mathscr{F}$}^{-1}\left\{\mathrm{rect}\left(\frac{(k_x-k_{x,0})}{D}\right)\tau_{\text{p}}(k_x,\omega)F(k_x,k_z,\omega)e^{-\sinh\left[\cosh^{-1}\left(\frac{c}{\omega n_b(\omega)}k_x\right)\right]z}\right\},
\end{align}
i.e., the inverse Fourier transform of a somewhat nasty function for the evanescent part where $\mathrm{rect}\left(\frac{(k_x-k_{x,0})}{D}\right)$ is the frequency window for the evanescent components. Evaluating the above inverse Fourier transform and computing the $\frac{1}{e}$ point of it would yield the depth (there is an integral of an exponentially decaying term here), I'm still working on this part...I'll update as it is completed.
