Waves in composites and metamaterials/Waves in layered media and point sources

The content of these notes is based on the lectures by Prof. Graeme W. Milton (University of Utah) given in a course on metamaterials in Spring 2007.

Previous Lecture

Recall from the previous lecture the governing equation for a TM wave.

$\text{(1)} \qquad \left[ \epsilon(z)~\cfrac{d }{d z}\left(\epsilon(z)^{-1}~\cfrac{d }{d z}\right) + k_z^2\right]~\varphi = 0$

where

$H_y(x, z) = \varphi(z)~e^{\pm i~k_x~x} = \tilde{H}_y(z)~e^{\pm i~k_x~x}$

and

$k_z^2 = \omega^2~\epsilon(z)~\mu(z) - k_x^2 ~.$

Also recall that we introduced a state vector $\mathbf{v}$ which is continuous across boundaries. The state vector is defined as

$\mathbf{v} := \begin{bmatrix} \varphi \\ \psi \end{bmatrix}$

where

$\psi := \cfrac{1}{i\omega\epsilon}~\cfrac{d \varphi}{d z} ~.$

Then,

$\text{(2)} \qquad \cfrac{d \mathbf{v}}{d z} = \begin{bmatrix} 0 & i\omega\epsilon \\ \cfrac{i~k_z^2}{\epsilon~\omega} & 0 \end{bmatrix}~\mathbf{v} = \overline{\mathbf{H}}~\mathbf{v} ~.$

The general solution of (2) is

$\text{(3)} \qquad \mathbf{v}(z) = A_+~e^{i~k_z~z}~\mathbf{n}_+ + A_-~e^{-i~k_z~z}~\mathbf{n}_-$

where $\mathbf{n}_+$ and $\mathbf{n}_-$ are the eigenvectors corresponding to the eigenvalues $i~k_z$ and $-i~k_z$, respectively.

We also found that the solution of the differential equation (2) can be expressed in terms of a propagator matrix $\overline{\mathbf{P}}$ such that

$\text{(4)} \qquad \mathbf{v}(z) = \overline{\mathbf{P}}(z, z')~\mathbf{v}(z')$

where

$\overline{\mathbf{P}}(z, z') := \tilde{\mathbf{n}} ~\overline{\mathbf{K}}(z - z')~\tilde{\mathbf{n}}^{-1} \equiv \tilde{\mathbf{n}}~e^{i\tilde{\mathbf{K}}(z - z')}~\tilde{\mathbf{n}}^{-1}$

where the columns of $\tilde{\mathbf{n}}$ are the eigenvectors of $\overline{\mathbf{H}}$ and

$\tilde{\mathbf{K}} = \begin{bmatrix} k_z & 0 \\ 0 & - k_z \end{bmatrix} ~;~~ e^{i\tilde{\mathbf{K}} z} = \begin{bmatrix} e^{i k_z z} & 0 \\ 0 & e^{-i k_z z} \end{bmatrix} ~.$

The matrix $\overline{\mathbf{P}}$ is called the propagator matrix or the transition matrix that related the fields at $z$ and $z'$.

Also, in a multilayered system (see Figure 1)

$\overline{\mathbf{P}}(-d_2, 0) = \overline{\mathbf{P}}(-d_2, -d_1) ~\overline{\mathbf{P}}(-d_1, 0)$

where $\overline{\mathbf{P}}(-d_2, -d_1)$ depends on $\epsilon_3, \mu_3$ and $\overline{\mathbf{P}}(-d_1, 0)$ depends on $\epsilon_2, \mu_2$.

Figure 1. Multilayered medium.

A special case

Let us now consider the particular case shown in Figure 2. The medium consists of three layers. Regions 1 and 3 are isotropic while the sandwiched region 2 is multilayered. The interface between regions 1 and 2 is located at $z = 0$ while the interface between regions 2 and 3 is located at $z = -d$. Let the propagator matrix of region 2 be $\overline{\mathbf{P}}(-d, 0)$. We want to find the reflection coefficient and the transmission coefficient of the system.

Figure 2. A multilayered medium sandwiched between two layers.

In reqion 1, the state vector is given by (see equation (3))

$\mathbf{v}_1(z) = A_{1_+}~e^{i~k_{z1}~z}~\mathbf{n}_{1_+} + A_{1_-}~e^{-i~k_{z1}~z}~\mathbf{n}_{1_-} ~.$

Define

$A_{1_+} := R~A_{1_-}$

where $R$ is like a scalar reflection coefficient. Then

$\text{(5)} \qquad \mathbf{v}_1(z) = R~A_{1_-}~e^{i~k_{z1}~z}~\mathbf{n}_{1_+} + A_{1_-}~e^{-i~k_{z1}~z}~\mathbf{n}_{1_-} ~.$

Proceeding as we did in the previous lecture, let us define matrices

$\tilde{\mathbf{n}}_1 := \begin{bmatrix} \mathbf{n}_{1_-} & \mathbf{n}_{1_+} \end{bmatrix} \qquad \text{and} \qquad e^{i\tilde{\mathbf{K}}_1 z} := \begin{bmatrix} e^{i k_{z1} z} & 0 \\ 0 & e^{-i k_{z1} z} \end{bmatrix} ~.$

Then equation (5) can be written as

$\text{(6)} \qquad { \mathbf{v}_1(z) = \tilde{\mathbf{n}}_1~e^{i\tilde{\mathbf{K}}_1 z}~ \begin{bmatrix} R \\ 1 \end{bmatrix}~A_{1_-} ~. }$

In region 3, there is only a transmitted wave. Therefore, the state vector is given by

$\text{(7)} \qquad \mathbf{v}_3(z) = A_{3_-}~e^{-i~k_{z3}~(z+d)}~\mathbf{n}_{3_-} ~.$

Define

$A_{3_-} := T~A_{1_-}$

where $T$ is a factor that acts like a scalar transmission coefficient. Then, equation (7) can be written in matrix form as

$\text{(8)} \qquad { \mathbf{v}_3(z) = \tilde{\mathbf{n}}_3~e^{i\tilde{\mathbf{K}}_3 (z+d)}~ \begin{bmatrix} 0 \\ T \end{bmatrix}~A_{1_-} }$

where

$\tilde{\mathbf{n}}_3 := \begin{bmatrix} \mathbf{n}_{3_-} & \mathbf{n}_{3_+} \end{bmatrix} \qquad \text{and} \qquad e^{i\tilde{\mathbf{K}}_3 z} := \begin{bmatrix} e^{i k_{z3} z} & 0 \\ 0 & e^{-i k_{z3} z} \end{bmatrix} ~.$

Since we have the propagator matrix for region 2, $\overline{\mathbf{P}}(-d,0)$, we can use it to connect regions 1 and 3. The continuity of the state vector across the interfaces implies that

$\text{(9)} \qquad \mathbf{v}_1(0) = \mathbf{v}_2(0) \quad \text{at}~z = 0\qquad \text{and} \qquad \mathbf{v}_2(-d) = \mathbf{v}_3(-d) \quad \text{at}~z = -d~.$

Also, using equation (4) we have

$\text{(10)} \qquad \mathbf{v}_2(-d) = \overline{\mathbf{P}}(-d, 0)~\mathbf{v}_2(0) ~.$

Therefore, using equations (9), we can write (10) as

$\text{(11)} \qquad { \mathbf{v}_3(-d) = \overline{\mathbf{P}}(-d, 0)~\mathbf{v}_1(0) ~. }$

From (6), at $z = 0$, we have

$\text{(12)} \qquad \mathbf{v}_1(0) = \tilde{\mathbf{n}}_1~\begin{bmatrix} R \\ 1 \end{bmatrix}~A_{1_-} ~.$

Also, from (8), at $z = -d$, we have

$\text{(13)} \qquad \mathbf{v}_3(-d) = \tilde{\mathbf{n}}_3~\begin{bmatrix} 0 \\ T \end{bmatrix}~A_{1_-} ~.$

Plugging into (11) gives

$\tilde{\mathbf{n}}_3~\begin{bmatrix} 0 \\ T \end{bmatrix}~A_{1_-} = \overline{\mathbf{P}}(-d, 0)~\tilde{\mathbf{n}}_1~\begin{bmatrix} R \\ 1 \end{bmatrix}~A_{1_-}$

or

$\text{(14)} \qquad { \begin{bmatrix} 0 \\ T \end{bmatrix} = \tilde{\mathbf{n}}^{-1}~\overline{\mathbf{P}}(-d, 0)~\tilde{\mathbf{n}}_1~ \begin{bmatrix} R \\ 1 \end{bmatrix} ~. }$

Equation (14) can then be solved to find the reflection and transmission coefficients $R$ and $T$.

Extension to waves in anisotropic layered media

In a layer medium where each of the layers is isotropic, the TE and TM waves are uncoupled at the interface. However, this is not true when each of the layers is anisotropic and we have to consider the full Maxwell's equations. The state vector approach can still be used for anisotropic media by choosing the variables such that they are continuous across interfaces.

Let us start with Maxwell's equations

$\boldsymbol{\nabla} \times \mathbf{E} = i \omega~\boldsymbol{\mu}\cdot\mathbf{H} ~;~~ \boldsymbol{\nabla} \times \mathbf{H} = -i \omega~\boldsymbol{\epsilon}\cdot\mathbf{E} ~.$

Recall that continuity of the fields requires that the tangential components of $\mathbf{E}$ and $\mathbf{H}$ be continuous across material interfaces. Therefore, an appropriate state vector for anisotropic media is

$\mathbf{v} := \begin{bmatrix} \mathbf{E}_s \\ \mathbf{H}_s \end{bmatrix}$

where $\mathbf{E}_s$ and $\mathbf{H}_s$ are the tangential components of $\mathbf{E}$ and $\mathbf{H}$ (i.e., the components on the surface normal to the $z$ direction).

Let us decompose the vector fields into a sum of the normal and tangential components:

$\mathbf{E} = \mathbf{E}_s + \mathbf{E}_z ~;~~ \mathbf{H} = \mathbf{H}_s + \mathbf{H}_z ~.$

The gradient operator can also be split along the same lines, i.e.,

$\boldsymbol{\nabla} = \boldsymbol{\nabla}_s + \frac{\partial }{\partial z}~\mathbf{e}_z ~;~~ \boldsymbol{\nabla}_s := \frac{\partial }{\partial x}~\mathbf{e}_x + \frac{\partial }{\partial y}~\mathbf{e}_y$

where $\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$ are the unit vectors in the $x, y, z$ directions, respectively. Let us express the tensors $\boldsymbol{\mu}$ and $\boldsymbol{\epsilon}$ in matrix form (with respect to the basis $\mathbf{e}_x, \mathbf{e}_y , \mathbf{e}_z$) as

$\boldsymbol{\mu} \equiv \begin{bmatrix} \boldsymbol{\mu}_{ss} & \boldsymbol{\mu}_{sz} \\ \boldsymbol{\mu}_{zs} & \mu_{zz} \end{bmatrix} ~;~~ \boldsymbol{\epsilon} \equiv \begin{bmatrix} \boldsymbol{\epsilon}_{ss} & \boldsymbol{\epsilon}_{sz} \\ \boldsymbol{\epsilon}_{zs} & \epsilon_{zz} \end{bmatrix}$

where $\boldsymbol{\mu}_{ss}, \boldsymbol{\epsilon}_{ss}$ are $2\times 2$ matrices, $\boldsymbol{\mu}_{sz}, \boldsymbol{\epsilon}_{sz}$ are $2\times 1$ matrices, $\boldsymbol{\mu}_{zs}, \boldsymbol{\epsilon}_{zs}$ are $1\times 2$ matrices, and $\mu_{zz}, \epsilon_{zz}$ are $1\times 1$ matrices, i.e., scalars.

Using the splits of the various quantities and the gradient operator, we can show ({\Red work this out}) that $\mathbf{E}_z, \mathbf{H}_z$ can be expressed in terms of $\mathbf{E}_s, \mathbf{H}_s$ as

$\mathbf{E}_z = -\cfrac{1}{i\omega \epsilon_{zz}}~\boldsymbol{\nabla}_s\times\mathbf{u}{\mathbf{H}_s} - \cfrac{1}{\epsilon_{zz}}~\boldsymbol{\epsilon}_{zs}~\mathbf{E}_s ~;~~ \mathbf{H}_z = \cfrac{1}{i\omega \mu_{zz}}~\boldsymbol{\nabla}_s\times\mathbf{u}{\mathbf{E}_s} - \cfrac{1}{\mu_{zz}}~\boldsymbol{\mu}{zs}~\mathbf{H}_s ~.$

After some further manipulations, the Maxwell equations may be expressed in matrix form as (see ~Chew95 for details)

$\text{(15)} \qquad { \frac{\partial }{\partial z} \begin{bmatrix} \mathbf{E}_s \\ \mathbf{H}_s \end{bmatrix} = \overline{\mathbf{H}} \begin{bmatrix} \mathbf{E}_s \\ \mathbf{H}_s \end{bmatrix} }$

where $\overline{\mathbf{H}}$ is a $4 \times 4$ matrix. If you compare equation (15) with (2), you will see that the form of the equations is the same, except that instead of a $2 \times 2$ matrix $\overline{\mathbf{H}}$ for the isotropic case, we now have a $4 \times 4$ matrix.

Plane wave expansion of sources in a homogeneous medium

It is often useful to expand sources in terms of plane waves so that the results of the previous lectures may be used directly on the basis of superposition. In this section we look at the expansion of point sources in terms of plane waves (for a homogeneous medium).

Let us look at the two-dimensional scalar wave equation (which can be used for acoustics, TE and TM waves, antiplane elasticity, etc.) In the presence of a point source, the wave equation has the form

$\text{(16)} \qquad \left[\frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} + k_\rho^2 \right]~\varphi(x, y) = -\delta(x)~\delta(y) ~.$

Assume that $k_\rho$ has a small positive imaginary part (it is a slightly lossy material), i.e.,

$k_\rho = k'_\rho + i k''_{\rho} ~.$

Expressed in cylindrical coordinates, equation (16) becomes (since the equation is symmetric about the origin)

$\text{(17)} \qquad \left[\frac{\partial^2 }{\partial \rho^2} + \cfrac{1}{\rho}~\frac{\partial^2 }{\partial \rho^2} + k_\rho^2 \right]~\varphi(\rho) = -\delta(\rho) ~;~~ \rho = \sqrt{x^2 + y^2} ~.$

The solution of (17) is

$\text{(18)} \qquad { \varphi(\rho) = \cfrac{i}{4}~H^{(1)}_0 (k_\rho~\rho) }$

where $H^{(1)}_0$ is a Hankel function of the first kind. \footnote{ Recall that a Hankel function of the first kind is defined as

$H^{(1)}_n(z) := J_n(z) + i~Y_n(z)$

where $J_n(z)$ is a Bessel function of the first kind and $Y_n(z)$ is a Bessel function of the second kind. }

We can also solve (16) using Fourier transforms. Let us assume that the function $\varphi(x, y)$ has a Fourier transform, i.e.,

$\text{(19)} \qquad \varphi(x,y) = \cfrac{1}{2\pi}\int_{-\infty}^{\infty} \widehat{\varphi}(k_x, y)~e^{i~k_x~x}~\text{d}k_x ~.$

Also note that

$\text{(20)} \qquad \delta(x) = \cfrac{1}{2\pi}\int_{-\infty}^{\infty} e^{i~k_x~x}~\text{d}k_x ~.$

Plugging (19) and (20) into (16) gives

$\cfrac{1}{2\pi}\int_{-\infty}^{\infty} \left[\frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} + k_\rho^2 \right]~ \left[\widehat{\varphi}(k_x, y) ~e^{i~k_x~x}\right]~\text{d}k_x = -\cfrac{1}{2\pi}\int_{-\infty}^{\infty} \delta(y)~e^{i~k_x~x}~\text{d}k_x$

or,

$\int_{-\infty}^{\infty} \left[-k_x^2 + \frac{\partial^2 }{\partial y^2} + k_\rho^2 \right]~ \widehat{\varphi}(k_x, y)~e^{i~k_x~x}~\text{d}k_x = -\int_{-\infty}^{\infty} \delta(y)~e^{i~k_x~x}~\text{d}k_x ~.$

Since the above equation holds for all values of $x$, the Fourier components must agree, i.e.,

$\left[\frac{\partial^2 }{\partial y^2} + k_\rho^2 - k_x^2 \right]~ \widehat{\varphi}(k_x, y) = -\delta(y) ~.$

Defining

$k_y^2 := k_\rho^2 - k_x^2$

we get

$\text{(21)} \qquad { \left[\frac{\partial^2 }{\partial y^2} + k_y^2 \right]~\widehat{\varphi}(k_x, y) = -\delta(y) ~. }$

Note that now the equation has a source only at $y = 0$. Away from the source (i.e., $y \ne 0$), the right hand side of (21) is zero, and the solution corresponds to the homogeneous part of the equation. Therefore,

$\text{(22)} \qquad \widehat{\varphi}(k_x, y) = A_1~e^{i~k_y~y} + A_2~e^{-i~k_y~y} \qquad y \ne 0 ~.$

This solution must be matched with the singularity at $y = 0$. This can be achieved by requiring that the solution have discontinuous second derivatives at $y = 0$. We then have ({\Red full explanation is needed here}), considering only waves that are damping away from the source rather that those growing exponentially,

$\text{(23)} \qquad \widehat{\varphi}(k_x, y) = \cfrac{i}{2~k_y}~e^{i~k_y~|y|} ~.$

Plugging (23) into (19) gives

$\text{(24)} \qquad { \varphi(x,y) = \cfrac{i}{4\pi}\int_{-\infty}^{\infty} \cfrac{1}{k_y}~ e^{i~k_x~x + i~k_y~|y|}~\text{d}k_x ~. }$

Equation (24) is a plane wave solution for the wave equation with a point source. So the point source has been converted into a sum of propagating plane waves and some evanescent terms.

Note that the denominator in (24) contains $k_y = \sqrt{k_\rho^2 - k_x^2}$. Hence, when $\pm k_\rho = k_x$ the solution blows up. Hence there are branch points at these locations as shown in Figure 3. In a lossless medium, $k''_\rho \rightarrow 0$ and these points appear as pole on the real $k_x$ axis. The integral in equation (24) can then be computed using the residual theorem. The region between the two poles is where waves are allowed to propagate in the $y$ direction while the region outside the poles is where these waves are evanescent.

Figure 3. Poles and integration path for plane wave solutions corresponding to a point source.

If we now compare the solutions (18) and (24), we have

$\text{(25)} \qquad H^{(1)}_0 (k_\rho~\rho) = \cfrac{1}{\pi}\int_{-\infty}^{\infty} \cfrac{1}{k_y}~ e^{i~k_x~x + i~k_y~|y|}~\text{d}k_x$

which provides a definition for the Hankel function.

Also, differentiating (16) with respect to $y$ and $x$, we get

$\left[\frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} + k_\rho^2 \right]~ \frac{\partial \varphi(x, y)}{\partial y} = -\delta(x)~\frac{\partial \delta(y)}{\partial y} \qquad \text{and} \qquad \left[\frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} + k_\rho^2 \right]~ \frac{\partial \varphi(x, y)}{\partial x} = -\delta(y)~\frac{\partial \delta(x)}{\partial x}$

Note that the products $\delta(x)~\frac{\partial \delta(y)}{\partial y}$ and $\delta(y)~\frac{\partial \delta(x)}{\partial x}$ correspond to dipole sources in the $y$ and $x$ directions, respectively.

Define

$\varphi_y(x, y) := \frac{\partial \varphi(x,y)}{\partial y} \qquad \text{and} \qquad \varphi_x(x, y) := \frac{\partial \varphi(x,y)}{\partial x} ~.$

Therefore, from (24), we have

${ \varphi_y(x,y) = -\cfrac{\text{sgn}(y)}{4\pi}\int_{-\infty}^{\infty} e^{i~k_x~x + i~k_y~|y|}~\text{d}k_x ~. \qquad \text{and} \qquad \varphi_x(x,y) = -\cfrac{1}{4\pi}\int_{-\infty}^{\infty} \cfrac{k_x}{k_y}~ e^{i~k_x~x + i~k_y~|y|}~\text{d}k_x ~. }$

These are the plane wave expansions of dipoles in the $y$ and $x$ directions respectively. Taking higher derivatives gives results from quadrupoles and other multipoles.

References

W. C. Chew. Waves and fields in inhomogeneous media. IEEE Press, New York, 1995.