Waves in composites and metamaterials/Continuum limit and propagator matrix

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 that we have been dealing with the TE equation

$\left[\mu (z)~{\cfrac {d}{dz}}\left(\mu (z)^{-1}~{\cfrac {d}{dz}}\right)+\omega ^{2}~\epsilon (z)~\mu (z)-k_{x}^{2}\right]~{\tilde {E}}_{y}=0$

where

$E_{y}(x,z)={\tilde {E}}_{y}(z)~e^{\pm i~k_{x}~x}~.$

For a a multilayered medium with $N$ layers, we found that in the $j$-th layer

${\text{(1)}}\qquad {\tilde {E}}_{yj}(z)=A_{j}~\left[\exp(-i~k_{zj}~z)+{\tilde {R}}_{j,j+1}~\exp(i~k_{zj}~z+2~i~k_{zj}~d_{j})\right]$

where ${\tilde {R}}_{j,j+1}$ is a generalized reflection coefficient. This coefficient can be obtained from a recursion relation of the form

${\text{(2)}}\qquad {\tilde {R}}_{j,j+1}={\cfrac {R_{j,j+1}+{\tilde {R}}_{j+1,j+2}~\exp[2~i~k_{z,j+1}~(d_{j+1}-d_{j})]}{1+R_{j,j+1}~{\tilde {R}}_{j+1,j+2}~\exp[2~i~k_{z,j+1}~(d_{j+1}-d_{j})]}}$

where

${\text{(3)}}\qquad R_{j,j+1}={\cfrac {\mu _{j+1}~k_{z,j}-\mu _{j}~k_{z,j+1}}{\mu _{j+1}~k_{z,j}+\mu _{j}~k_{z,j+1}}}$

is the Fresnel reflection coefficient for TE waves. Equation (3) may also be written as

${\text{(4)}}\qquad R_{j,j+1}={\cfrac {{\cfrac {k_{z,j}}{\mu _{j}}}-{\cfrac {k_{z,j+1}}{\mu _{j+1}}}}{{\cfrac {k_{z,j}}{\mu _{j}}}+{\cfrac {k_{z,j+1}}{\mu _{j+1}}}}}~.$

We will now proceed to determine the generalized reflection coefficient in the continuum limit.

Continuum Limit

Consider a medium where the layer thickness if $\Delta$. Denote the reflection coefficient and the generalized reflection coefficient at the interface $z=-d_{j}$ as $R(z)$ and ${\tilde {R}}(z)$, respectively. Also denote the phase velocity $k_{z,j+1}$ just below the interface as $k_{z}(z-\Delta /2)$ and the permeability $\mu _{j+1}$ as $\mu (z-\Delta /2)$. \footnote {This implies that we are measuring the phase velocity and the permeability at the center of the layer. However, this is not strictly necessary and we could alternatively measure these quantities at $z-\Delta$.}

Then, equation (2) can be written as

${\text{(5)}}\qquad {\tilde {R}}(z)={\cfrac {R(z)+{\tilde {R}}(z-\Delta )~\exp \left[2~i~k_{z}\left(z-\Delta /2\right)~\Delta \right]}{1+R(z)~{\tilde {R}}(z-\Delta )~\exp \left[2~i~k_{z}\left(z-\Delta /2\right)~\Delta \right]}}$

where

$R(z)={\cfrac {{\cfrac {k_{z}(z+\Delta /2)}{\mu (z+\Delta /2)}}-{\cfrac {k_{z}(z-\Delta /2)}{\mu (z-\Delta /2)}}}{{\cfrac {k_{z}(z+\Delta /2)}{\mu (z+\Delta /2)}}+{\cfrac {k_{z}(z-\Delta /2)}{\mu (z-\Delta /2)}}}}={\cfrac {{\tilde {k}}(z+\Delta /2)-{\tilde {k}}(z-\Delta /2)}{{\tilde {k}}(z+\Delta /2)+{\tilde {k}}(z-\Delta /2)}}$

with

${\tilde {k}}:={\cfrac {k_{z}}{\mu }}~.$

Expanding in Taylor series about $z$ and ignoring higher order terms, we get

${\text{(6)}}\qquad R(z)\approx {\cfrac {{\tilde {k}}(z)+\Delta /2~{\tilde {k}}^{'}(z)-{\tilde {k}}(z)+\Delta /2~{\tilde {k}}^{'}(z)}{{\tilde {k}}(z)+\Delta /2~{\tilde {k}}^{'}(z)+{\tilde {k}}(z)-\Delta /2~{\tilde {k}}^{'}(z)}}={\cfrac {\Delta ~{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}~.$

Similarly, igoring powers $\Delta ^{2}$ and higher, we get

${\text{(7)}}\qquad \exp \left[2~i~k_{z}\left(z-\Delta /2\right)~\Delta \right]\approx 1+2~\Delta ~i~k_{z}(z)$

and

${\text{(8)}}\qquad {\tilde {R}}(z-\Delta )\approx {\tilde {R}}(z)-\Delta ~{\tilde {R}}^{'}(z)~.$

Plugging the expansions (7) and (8) into (5) gives

${\text{(9)}}\qquad {\tilde {R}}(z)\approx {\cfrac {R(z)+\left[{\tilde {R}}(z)-\Delta ~{\tilde {R}}^{'}(z)\right]~\left[1+2~\Delta ~i~k_{z}(z)\right]}{1+R(z)~\left[{\tilde {R}}(z)-\Delta ~{\tilde {R}}^{'}(z)\right]~\left[1+2~\Delta ~i~k_{z}(z)\right]}}~.$

Substituting (6) into (9) and dropping terms containing $\Delta ^{2}$ and higher gives

${\text{(10)}}\qquad {\tilde {R}}(z)\approx {\cfrac {{\cfrac {\Delta ~{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}+{\tilde {R}}(z)-\Delta ~{\tilde {R}}^{'}(z)+2~\Delta ~i~k_{z}(z)~{\tilde {R}}(z)}{1+{\cfrac {\Delta ~{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}~{\tilde {R}}(z)}}~.$

If we assume that $\Delta ~{\tilde {k}}^{'}$ is small such that the denominator can be expanded in a series, we get

${\text{(11)}}\qquad {\tilde {R}}(z)\approx \left[{\cfrac {\Delta ~{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}+{\tilde {R}}(z)-\Delta ~{\tilde {R}}^{'}(z)+2~\Delta ~i~k_{z}(z)~{\tilde {R}}(z)\right]\left[1-{\cfrac {\Delta ~{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}~{\tilde {R}}(z)+O(\Delta ^{2})\right]~.$

After expanding and ingoring terms containing $\Delta ^{2}$, we get

${\text{(12)}}\qquad {\tilde {R}}(z)\approx {\cfrac {\Delta ~{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}+{\tilde {R}}(z)-\Delta ~{\tilde {R}}^{'}(z)+2~\Delta ~i~k_{z}(z)~{\tilde {R}}(z)-{\cfrac {\Delta ~{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}~[{\tilde {R}}(z)]^{2}$

or,

${\text{(13)}}\qquad {{\tilde {R}}^{'}(z)=2~i~k_{z}(z)~{\tilde {R}}(z)+{\cfrac {{\tilde {k}}^{'}(z)}{2~{\tilde {k}}(z)}}\left\{1-[{\tilde {R}}(z)]^{2}\right\}~.}\qquad {\text{Ricotti equation}}$

We thus get an equation that gives a continuous representation of the generalized reflection coefficient ${\tilde {R}}(z)$. Equation (13) can be solved numerically using the Runge-Kutta method.

For example, in the stuation shown in Figure 1, the generalized reflectivity coefficient at the point $z_{1}$ is ${\tilde {R}}(z_{1})$ while that at point $z_{0}$ is 0. If we wish to determine the value of ${\tilde {R}}(z_{i})$ at a point inside the smoothly varying layer, then one possibility is to assume that $\mu (z)$ and $\epsilon (z)$ is constant for $z>z_{i}$ and compute the value of ${\tilde {R}}$ in the usual manner.

Figure 1. Reflectivity in a smoothly graded layered material.

There can also be a situation where there are a few isolated strong discontinuities inside the graded layer as shown in Figure 2. If there is a discontinuity at $z_{c}$, we can use the discrete solution with layer thickness 0 at the discontinuity.

Figure 2. Reflectivity in a smoothly graded material with a strong discontinuity.

Then, from (2), at the discontinuity

${\text{(14)}}\qquad {\tilde {R}}(z_{c}^{+})={\cfrac {R(z_{c})+{\tilde {R}}(z_{c}^{-})}{1+R(z_{c})~{\tilde {R}}(z_{c}^{-})}}~.$

Also, from (4)

${\text{(15)}}\qquad R(z_{c})={\cfrac {{\cfrac {k_{z}^{+}}{\mu ^{+}}}-{\cfrac {k_{z}^{-}}{\mu ^{-}}}}{{\cfrac {k_{z}^{+}}{\mu ^{+}}}+{\cfrac {k_{z}^{-}}{\mu ^{-}}}}}~.$

Hence we can find the generalized reflection coefficients at isolated discontinuities within the material.

Determining the coefficients $A_{j}$

Recall equation (1):

${\tilde {E}}_{yj}(z)=A_{j}~\left[\exp(-i~k_{zj}~z)+{\tilde {R}}_{j,j+1}~\exp(i~k_{zj}~z+2~i~k_{zj}~d_{j})\right]~.$

So far we have determine the value of ${\tilde {R}}$ is this equation. But how do we determine the coefficients $A_{j}$ in multilayered media?

Let us start with the coefficients for a single layer that we determined in the previous lecture. We had

${\text{(16)}}\qquad A_{2}={\cfrac {T_{12}~A_{1}~\exp[i~(k_{z1}-k_{z2})~d_{1}]}{1-R_{21}~R_{23}~\exp[2~i~k_{z2}~(d_{2}-d_{1})]}}$

where

$T_{12}={\cfrac {2~\mu _{2}~k_{z1}}{\mu _{2}~k_{z1}+\mu _{1}~k_{z2}}}~.$

We can rewrite (16) as

${\text{(17)}}\qquad A_{2}~\exp(i~k_{z2}~d1)={\cfrac {T_{12}~A_{1}~\exp(i~k_{z1}~d_{1})}{1-R_{21}~R_{23}~\exp[2~i~k_{z2}~(d_{2}-d_{1})]}}$

Using the same arguments as before, we can generalize (17) to a medium with $N$ layers. Thus, for the $j$-th layer, we have

${\text{(18)}}\qquad A_{j}~\exp(i~k_{z,j}~d_{j-1})={\cfrac {T_{j-1,j}~A_{j-1}~\exp(i~k_{z,j-1}~d_{j-1})}{1-R_{j,j-1}~{\tilde {R}}_{j,j+1}~\exp[2~i~k_{z,j}~(d_{j}-d_{j-1})]}}~.$

Define

$S_{j-1,j}:={\cfrac {T_{j-1,j}}{1-R_{j,j-1}~{\tilde {R}}_{j,j+1}~\exp[2~i~k_{z,j}~(d_{j}-d_{j-1})]}}~.$

Then we can write (18) as

${\text{(19)}}\qquad {\begin{aligned}A_{j}~\exp(i~k_{z,j}~d_{j-1})&=S_{j-1,j}~A_{j-1}~\exp(i~k_{z,j-1}~d_{j-1})\\&=S_{j-1,j}~A_{j-1}~\exp(i~k_{z,j-1}~d_{j-2})~\exp[i~k_{z,j-1}~(d_{j-1}-d_{j-2})]~.\end{aligned}}$

The second of equations (19) gives us a recurrence relation that can be used to compute the other $A_{j}$ s. Thus, we can write

${\text{(20)}}\qquad {\begin{aligned}A_{j}~\exp(i~k_{z,j}~d_{j-1})&=A_{j-1}~\exp(i~k_{z,j-1}~d_{j-2})~S_{j-1,j}~\exp[i~k_{z,j-1}~(d_{j-1}-d_{j-2})]\\&=A_{j-2}~\exp(i~k_{z,j-2}~d_{j-3})~S_{j-2,j-1}~\exp[i~k_{z,j-2}~(d_{j-2}-d_{j-3})]~S_{j-1,j}~\exp[i~k_{z,j-1}~(d_{j-1}-d_{j-2})]\\&=\dots \\&=A_{1}~\exp(i~k_{z1}~d_{1})~\prod _{m=1}^{j-1}S_{m,m+1}~\exp[i~k_{z,m}~(d_{m}-d_{m-1})]~.\end{aligned}}$

We can introduce a generalized transmission coefficient

${\tilde {T}}_{1N}:=\prod _{m=1}^{N-1}S_{m,m+1}~\exp[i~k_{z,m}~(d_{m}-d_{m-1})]~.$

Then,

${A_{N}~\exp(i~k_{z,N}~d_{N-1})={\tilde {T}}_{1N}~A_{1}~\exp(i~k_{z1}~d_{1})~.}$

So the downgoing wave amplitude in region $N$ at $z=-d_{n+1}$ is ${\tilde {T}}_{1N}$ times the downgoing amplitude in region 1 ($z=-d_{1}$). Due to the products involved in the above relation, a continuum extension of this formula is not straightforward.

State equations and Propagator matrix

The propagator matrix relates the fields at two points in a multilayered medium. This matrix is also known as the transition matrix or the transfer matrix.

Let us examine the propagation matrix for a TM wave. Recall the governing equation for a TM wave:

${\text{(21)}}\qquad \left[\epsilon (z)~{\cfrac {d}{dz}}\left(\epsilon (z)^{-1}~{\cfrac {d}{dz}}\right)+\omega ^{2}~\epsilon (z)~\mu (z)-k_{x}^{2}\right]~{\tilde {H}}_{y}=0$

where

$H_{y}(x,z)={\tilde {H}}_{y}(z)~e^{\pm i~k_{x}~x}~.$

Define $\varphi :={\tilde {H}}_{y}$. Also,

$k_{z}^{2}=\omega ^{2}~\epsilon (z)~\mu (z)-k_{x}^{2}~.$

Therefore, (21) can be written as

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

To reduce (22) to a first order differential equation, introduce the quantity

${\text{(23)}}\qquad \psi :={\cfrac {1}{i\omega \epsilon }}~{\cfrac {d\varphi }{dz}}\qquad \implies \qquad {\cfrac {d\varphi }{dz}}=i\omega \epsilon ~\psi ~.$

Clearly, $\psi$ has to be continuous across the interface for the differential equation (22) to be satisfied.

Plugging (23) into (22) gives

${\text{(24)}}\qquad {\cfrac {d\psi }{dz}}={\cfrac {i~k_{z}^{2}}{\epsilon ~\omega }}~\varphi ~.$

Therefore, (23)$_{2}$ and (24) for a system of differential equations which can be written as

${\text{(25)}}\qquad {\cfrac {d}{dz}}{\begin{bmatrix}\varphi \\\psi \end{bmatrix}}={\begin{bmatrix}0&i\omega \epsilon \\{\cfrac {i~k_{z}^{2}}{\epsilon ~\omega }}&0\end{bmatrix}}~{\begin{bmatrix}\varphi \\\psi \end{bmatrix}}~.$

Define

$\mathbf {v} :={\begin{bmatrix}\varphi \\\psi \end{bmatrix}}~;~~{\overline {\mathbf {H} }}:={\begin{bmatrix}0&i\omega \epsilon \\{\cfrac {i~k_{z}^{2}}{\epsilon ~\omega }}&0\end{bmatrix}}~.$

Then, equations (25) can be written as

${\text{(26)}}\qquad {\cfrac {d\mathbf {v} }{dz}}={\overline {\mathbf {H} }}~\mathbf {v} ~.$

If ${\overline {\mathbf {H} }}$ is constant, particular solutions to (26) can be sought of the form

${\text{(27)}}\qquad \mathbf {v} =\mathbf {v} _{0}~e^{\lambda ~z}~.$

Plugging (27) into (25) leads to the eigenvalue problem

$({\overline {\mathbf {H} }}-\lambda ~\mathbf {1} )~\mathbf {v} _{0}=\mathbf {0} ~.$

Solutions exist only if

$\det({\overline {\mathbf {H} }}-\lambda ~\mathbf {1} )=0\qquad \implies \qquad \lambda ^{2}=-k_{z}^{2}\qquad \implies \qquad \lambda =\pm i~k_{z}~.$

Therefore, the general solution of (26) is

${\text{(28)}}\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.

Equation (28) can be written more compactly in the form

${\text{(29)}}\qquad \mathbf {v} (z)={\begin{bmatrix}\mathbf {n} _{+}&\mathbf {n} _{-}\end{bmatrix}}{\begin{bmatrix}e^{i~k_{z}~z}&0\\0&e^{-i~k_{z}~z}\end{bmatrix}}{\begin{bmatrix}A_{+}\\A_{-}\end{bmatrix}}$

or,

${\text{(30)}}\qquad \mathbf {v} (z)={\tilde {\mathbf {n} }}~{\overline {\mathbf {K} }}(z)~{\tilde {\mathbf {a} }}$

where

${\tilde {\mathbf {n} }}:={\begin{bmatrix}\mathbf {n} _{+}&\mathbf {n} _{-}\end{bmatrix}}~;~~{\overline {\mathbf {K} }}(z):={\begin{bmatrix}e^{i~k_{z}~z}&0\\0&e^{-i~k_{z}~z}\end{bmatrix}}~;~~{\tilde {\mathbf {a} }}:={\begin{bmatrix}A_{+}\\A_{-}\end{bmatrix}}~.$

Note that, for a point $z'$ that is different from $z$,

${\overline {\mathbf {K} }}(z-z')={\begin{bmatrix}e^{i~k_{z}~(z-z')}&0\\0&e^{-i~k_{z}~(z-z')}\end{bmatrix}}={\begin{bmatrix}e^{i~k_{z}~z}&0\\0&e^{-i~k_{z}~z}\end{bmatrix}}{\begin{bmatrix}e^{-i~k_{z}~z'}&0\\0&e^{i~k_{z}~z'}\end{bmatrix}}={\overline {\mathbf {K} }}(z)~{\overline {\mathbf {K} }}(-z')~.$

Also,

${\overline {\mathbf {K} }}(z-z')~{\overline {\mathbf {K} }}(z')={\begin{bmatrix}e^{i~k_{z}~(z-z')}&0\\0&e^{-i~k_{z}~(z-z')}\end{bmatrix}}{\begin{bmatrix}e^{-i~k_{z}~z'}&0\\0&e^{i~k_{z}~z'}\end{bmatrix}}={\overline {\mathbf {K} }}(z)~.$

Therefore we can write (30) in the form

$\mathbf {v} (z)={\tilde {\mathbf {n} }}~{\overline {\mathbf {K} }}(z-z')~{\overline {\mathbf {K} }}(z')~{\tilde {\mathbf {a} }}={\tilde {\mathbf {n} }}~{\overline {\mathbf {K} }}(z-z')~{\tilde {\mathbf {n} }}^{-1}~~{\tilde {\mathbf {n} }}~{\overline {\mathbf {K} }}(z')~{\tilde {\mathbf {a} }}$

or,

${\text{(31)}}\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}~.}$

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

In a multilayered system (see Figure~3), since the vector $\mathbf {v}$ is discontinuous, we can show that

${\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 3. Multilayered medium.

References

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