Waves in composites and metamaterials/Airy solution and WKB solution

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.

Introduction

Recall from the previous lecture that we assumed that the permittivity and permeability are scalars and are locally isotropic though not globally so. Then we may write

${\boldsymbol {\epsilon }}=\epsilon (x_{3})~{\boldsymbol {\mathit {1}}}\quad {\text{and}}\quad {\boldsymbol {\mu }}=\mu (x_{3})~{\boldsymbol {\mathit {1}}}~.$

The TE (transverse electric field) equations are given by

${\text{(1)}}\qquad {\overline {\boldsymbol {\nabla }}}\cdot \left({\cfrac {1}{\mu }}~{\overline {\boldsymbol {\nabla }}}E_{1}\right)+\omega ^{2}~\epsilon ~E_{1}=0$

where ${\overline {\boldsymbol {\nabla }}}$ represents the two-dimensional gradient operator. Equation (1) can also be written as

${\text{(2)}}\qquad \left[{\frac {\partial }{\partial x_{2}^{2}}}+\mu (x_{3})~{\frac {\partial }{\partial x_{3}}}\left({\cfrac {1}{\mu (x_{3})}}~{\frac {\partial }{\partial x_{3}}}\right)+\omega ^{2}~\epsilon (x_{3})~\mu (x_{3})\right]~E_{1}=0$

which admits solutions of the form

$E_{1}(x_{2},x_{3})={\tilde {E}}_{1}(x_{3})~e^{\pm i~k_{2}~x_{2}}$

and equation (2) then becomes an ODE:

${\text{(3)}}\qquad \left[\mu (x_{3})~{\cfrac {d}{dx_{3}}}\left({\cfrac {1}{\mu (x_{3})}}~{\cfrac {d}{dx_{3}}}\right)+\omega ^{2}~\epsilon (x_{3})~\mu (x_{3})-k_{2}^{2}\right]~{\tilde {E}}_{1}=0~.$

The quantity

$k_{3}^{2}:=\omega ^{2}~\epsilon (x_{3})~\mu (x_{3})-k_{2}^{2}$

can be less than zero, implying that $k_{3}$ may be complex. Also, at the boundary, both ${\tilde {E}}_{1}$ and $1/\mu \partial {\tilde {E}}_{1}/\partial x_{3}$ must be continuous.

TE waves in a non-magnetic medium

For a non-magnetic medium, $\mu$ is constant and we can write (3) as

${\text{(4)}}\qquad \left[{\cfrac {d^{2}}{dx^{2}}}+\omega ^{2}~\epsilon ~\mu -k_{2}^{2}\right]~{\tilde {E}}_{1}=0\qquad {\text{where}}\quad x\equiv x_{3}~.$

Permittivity varies linearly with x

If the permittivity varies linearly with $x$, then we may write

$\epsilon (x)=a+b~x$

where $a$ and $b$ are constants. Plugging this into (4) we get

${\text{(5)}}\qquad \left[{\cfrac {d^{2}}{dx^{2}}}+A+B~x\right]~{\tilde {E}}_{1}=0\qquad {\text{where}}\quad A:=\omega ^{2}~\mu ~a-k_{2}^{2}~;~~B:=\omega ^{2}~\mu ~b~.$

Let us assume that $B>0$ (this is not strictly necessary, but simplifies things for our present analysis). Let us introduce a change of variables

$\eta =B^{1/3}\left(x+{\cfrac {A}{B}}\right)~.$

Then (5) becomes

${\text{(6)}}\qquad \left[{\cfrac {d^{2}}{d\eta ^{2}}}+\eta \right]~{\tilde {E}}_{1}=0~.$

Equation (6) is called the Airy equation. The solution of this equation is

${\tilde {E}}_{1}(\eta )=C_{1}~{\text{Ai}}(-\eta )+C_{2}~{\text{Bi}}(-\eta )$

where ${\text{Ai}}$ and ${\text{Bi}}$ are Airy functions of the first and second kind (see Abram72 for details.) A plot of the behavior of the two Airy functions as a function of real $-\eta$ is shown in Figure~1.

Figure 1. The Airy functions ${\text{Ai}}$ and ${\text{Bi}}$ as functions of $-\eta$.

As $x\rightarrow -\infty$ (i.e., as $-\eta \rightarrow \infty$), the Airy functions asymptotically approach the values

${\begin{aligned}{\text{Ai}}(-\eta )&\sim {\cfrac {1}{2}}~\pi ^{-1/2}~(-\eta )^{-1/4}~e^{-2/3~(-\eta )^{3/2}}\\{\text{Bi}}(-\eta )&\sim \pi ^{-1/2}~(-\eta )^{-1/4}~e^{2/3~(-\eta )^{3/2}}~.\end{aligned}}$

So ${\text{Ai}}(-\eta )$ corresponds to an exponentially decaying wave as $|x|\rightarrow \infty$ and ${\text{Bi}}(-\eta )$ corresponds to an exponentially increasing waves at $|x|\rightarrow \infty$. A schematic of the situation is shown in Figure 2.

Figure 2. The region where the TE wave is exponentially damped.

If there are no sources in the region $x<0$ then the solution ${\text{Bi}}(-\eta )$ is unphysical which implies that $C_{2}=0$. Therefore,

${\text{(7)}}\qquad {\tilde {E}}_{1}(\eta )=C_{1}~{\text{Ai}}(-\eta )~.$

Now, as $x\rightarrow \infty$ (i.e., as $-\eta \rightarrow -\infty$), the Airy function ${\text{Ai}}(-\eta )$ takes the asymptotic form

${\text{(8)}}\qquad {\text{Ai}}(-\eta )\sim \pi ^{-1/2}~\eta ^{-1/4}~\sin \left({\cfrac {2}{3}}~\eta ^{3/2}+{\cfrac {\pi }{4}}\right)~.$

This is a superposition of right and left travelling waves (because the sine can be decomposed into two exponentials one of which corresponds to a wave travelling in one direction and the seconds to a wave travelling in the opposite direction.)

The Wentzel-Kramers-Brillouin (WKB) method

If we don't assume any particular linear variation of the permittivity $\epsilon (x)$, we can use the WKB method to arrive at a solution for high frequency waves.

The WKB method is a high frequency method for obtaining solutions to one-dimensional (time-independent) wave equations of the form

${\text{(9)}}\qquad {\cfrac {d^{2}\varphi }{dx^{2}}}+k_{3}^{2}(x)~\varphi (x)=0~.$

Recall from (1) that the TE equation in a nonmagnetic medium is

${\cfrac {d^{2}{\tilde {E}}_{1}}{dx^{2}}}+(\omega ^{2}~\epsilon ~\mu -k_{2}^{2})~{\tilde {E}}_{1}=0~.$

Clearly this equation can be written in form (9) by setting

$\varphi ={\tilde {E}}_{1}~;~~k_{3}^{2}(x)=\omega ^{2}~\epsilon ~\mu -k_{2}^{2}~.$

Recall also that the TM equation is

${\text{(10)}}\qquad \epsilon ~{\cfrac {d}{dx}}\left(\epsilon ^{-1}~{\cfrac {d{\tilde {H}}_{1}}{dx}}\right)++k_{3}^{2}(x)~{\tilde {H}}_{1}=0~.$

Equation (11) can also be reduced to the form (9). The procedure is as follows. Let us first set $\psi ={\tilde {H}}_{1}$ to get

${\text{(11)}}\qquad \epsilon ~{\cfrac {d}{dx}}\left(\epsilon ^{-1}~{\cfrac {d\psi }{dx}}\right)+k_{3}^{2}(x)~\psi =0~.$

After expanding (11) we get

${\text{(12)}}\qquad {\cfrac {d^{2}\psi }{dx^{2}}}=\epsilon ^{-1}~{\cfrac {d\epsilon }{dx}}~{\cfrac {d\psi }{dx}}-k_{3}^{2}~\psi =0~.$

Define

${\text{(13)}}\qquad \varphi (x):=\epsilon ^{-1/2}(x)~\psi (x)~.$

Differentiating (13) twice, we get

${\text{(14)}}\qquad {\cfrac {d^{2}\varphi }{dx^{2}}}={\cfrac {d^{2}}{dx^{2}}}\left(\epsilon ^{-1/2}\right)~\psi +\epsilon ^{-1/2}~{\cfrac {d^{2}\psi }{dx^{2}}}+2~{\cfrac {d}{dx}}\left(\epsilon ^{-1/2}\right)~{\cfrac {d\psi }{dx}}~.$

Substituting (12), (13) into (14) we have

${\text{(15)}}\qquad {\cfrac {d^{2}\varphi }{dx^{2}}}={\cfrac {d^{2}}{dx^{2}}}\left(\epsilon ^{-1/2}\right)~(\epsilon ^{1/2}~\phi )+\epsilon ^{-3/2}~{\cfrac {d\epsilon }{dx}}~{\cfrac {d\psi }{dx}}-k_{3}^{2}~\epsilon ^{-1/2}~\psi -\epsilon ^{-3/2}~{\cfrac {d\epsilon }{dx}}~{\cfrac {d\psi }{dx}}$

or,

${\text{(16)}}\qquad {\cfrac {d^{2}\varphi }{dx^{2}}}+\left[k_{3}^{2}-\epsilon ^{1/2}~{\cfrac {d^{2}}{dx^{2}}}\left(\epsilon ^{-1/2}\right)\right]~\varphi =0~.$

Equation (16) has the same form as (9).

At this stage recall that

$k_{3}^{2}=\omega ^{2}~\epsilon ~\mu -k_{2}^{2}~.$

Let us assume that $k_{2}$ is proportional to $\omega$ which implies that $k_{3}$ is also proportional to $\omega$, i.e.,

${\text{(17)}}\qquad k_{3}^{2}(x)=\omega ^{2}~s^{2}(x)$

where $s(x)$ is independent of $\omega$.

In equation (16), if $\omega$ is large, then $k_{3}$ will dominate and we will end up with exactly the same equation as (9), provided variations in $\epsilon$ are smooth (and we don't get large jumps in its derivatives).

Let us now try to solve (9). When $k_{3}$ is constant, the solution of the equation is a traveling wave. If we assume that $k_{3}$ varies slowly with $x$, we can try to get solutions of the form

${\text{(18)}}\qquad \varphi (x)=A~e^{i\omega ~\tau (x)}$

and examine the phase $\tau (x)$ rather than the solution $\varphi (x)$. Differentiating (18), we get

${\text{(19)}}\qquad \varphi '(x)=i\omega ~\tau '(x)~A~e^{i\omega ~\tau (x)}~;~~\varphi ''(x)=\left[i\omega ~\tau ''(x)-\omega ^{2}~(\tau '(x))^{2}\right]~A~e^{i\omega ~\tau (x)}~.$

Plugging (19) into (9), we get

${\text{(20)}}\qquad i\omega ~\tau ''(x)-\omega ^{2}~(\tau '(x))^{2}+k_{3}^{2}(x)=0~.$

If we assume that $k_{3}^{2}>0$ (i.e., $k_{3}$ is real) we can simplify the analysis slightly at this stage (even though this is not strictly necessary).

For large $\omega$, i.e., $\omega \gg 1$, we can seek a perturbation solution of the form

${\text{(21)}}\qquad \tau (x)=\tau _{0}(x)+{\cfrac {1}{\omega }}~\tau _{1}(x)+{\cfrac {1}{\omega ^{2}}}~\tau _{2}(x)+\dots ~.$

Plugging (21) into (20) and using (17) we get

${\text{(22)}}\qquad \left[i\omega ~\tau _{0}''(x)+i~\tau _{1}''(x)+{\cfrac {i}{\omega }}~\tau _{2}''(x)+\dots \right]-\left[\omega ~\tau _{0}'(x)+\tau _{1}'(x)+{\cfrac {1}{\omega }}~\tau _{2}'(x)+\dots \right]^{2}+\omega ^{2}~s^{2}(x)=0$

or,

${\text{(23)}}\qquad \left[{\cfrac {i}{\omega }}~\tau _{0}''(x)+{\cfrac {i}{\omega ^{2}}}~\tau _{1}''(x)+{\cfrac {i}{\omega ^{3}}}~\tau _{2}''(x)+\dots \right]-\left[\tau _{0}'(x)+{\cfrac {1}{\omega }}~\tau _{1}'(x)+{\cfrac {1}{\omega ^{2}}}~\tau _{2}'(x)+\dots \right]^{2}+s^{2}(x)=0~.$

For large $\omega$, equation (23) reduces to

${\text{(24)}}\qquad -[\tau _{0}'(x)]^{2}+s^{2}(x)=0\qquad {\text{or}}\qquad [\tau _{0}'(x)]^{2}=s^{2}(x)\qquad \qquad {\text{(Eikonal equation)}}~.$

Therefore,

${\text{(25)}}\qquad \tau _{0}'(x)=\pm s(x)~.$

Integrating (25) from an arbitrary point $x_{0}$ to $x$, we get

${\text{(26)}}\qquad \tau _{0}(x)=\pm \int _{x_{0}}^{x}s(y)~{\text{d}}y+C_{0_{\pm }}.$

where $C_{0_{\pm }}$ depends on the sign of the integral.

Next, collecting terms of order $\omega$ in equation (22), we get

${\text{(27)}}\qquad i\omega ~\tau _{0}''(x)-2~\omega ~\tau _{0}'(x)~\tau _{1}'(x)=0~.$

Substituting (25) into (27) we get

$\pm i\omega ~s'(x)\mp 2~\omega ~s(x)~\tau _{1}'(x)=0$

or,

${\text{(28)}}\qquad {\cfrac {i}{2}}~{\cfrac {s'(x)}{s(x)}}=\tau _{1}'(x)~.$

Integrating (28) we get

${\text{(29)}}\qquad \tau _{1}(x)={\cfrac {i}{2}}~\ln[s(x)]+C_{1}=i\ln[{\sqrt {s(x)}}]+C_{1}~.$

Plugging (26) and (29) into (21) (and ignoring terms containing powers of $\omega ^{2}$ and higher) we get

${\text{(30)}}\qquad \tau (x)=\pm \int _{x_{0}}^{x}s(y)~{\text{d}}y+{\cfrac {i}{\omega }}~\ln[{\sqrt {s(x)}}]+C_{\pm }~.$

This implies that the solution (18) has the form

${\text{(31)}}\qquad {\varphi (x)={\cfrac {A_{+}}{\sqrt {s(x)}}}~\exp \left(i\omega ~\int _{x_{0}}^{x}s(y)~{\text{d}}y\right)+{\cfrac {A_{-}}{\sqrt {s(x)}}}~\exp \left(-i\omega ~\int _{x_{0}}^{x}s(y)~{\text{d}}y\right)~.}$

Equation (31) is the WKB solution assuming $k_{3}^{2}>0$. Note that when $s(x)=0$, a solution does not exist.

Also note that since $k_{3}^{2}$ is proportional to $\omega ^{2}$,

${\text{(32)}}\qquad |k_{3}^{2}|\gg |i\omega ~\tau _{0}''|\qquad {\text{for large}}~\omega ~.$

Therefore,

$\omega ^{2}~s^{2}(x)\gg \omega ~s'(x)\qquad \implies \qquad \omega ~s(x)\gg {\cfrac {\omega ~s'(x)}{\omega ~s(x)}}={\cfrac {d}{dx}}\left[\ln(\omega ~s(x))\right]$

or,

$k_{3}(x)\gg {\cfrac {d}{dx}}\left[\ln(k_{3}(x))\right]~.$

Therefore, the restriction is that $\omega$ is large and that $k_{3}$ is smooth with respect to $x$.

Now, consider for example the profile shown in Figure 3. In region I, the WKB solution is valid since $k_{3}^{2}>0$. At the point where the profile meets the $x$ axis, a solution does not exist since $s(x)=0$. However, if the profile is smooth enough, we can assume that $k_{3}(x)$ is linear and we can use the Airy solution for the region II around this point. When the profile goes below the $x$ axis, $k_{3}^{2}<0$. However, the WKB solution is valid in this region (III) too as equation 32 can still be satisfied with $s(x)=i~\alpha (x)$.

Figure 3. Regions of validity of the linear solution and the WKB solution for large $\omega$.

There is an area of overlap between the regions where the WKB solution is valid and the region where the Airy solution is valid. In fact, the unknown parameters in the two solutions can be determined by matching the solutions at points in this region of overlap.

To do this, let $\zeta$ be the point on the $x$-axis where $s(x)=0$. Then, in region I, the solution is

${\text{(33)}}\qquad \varphi _{I}(x)={\cfrac {A_{+}}{\sqrt {s(x)}}}~\exp \left(i\omega ~\int _{\zeta }^{x}s(y)~{\text{d}}y\right)+{\cfrac {A_{-}}{\sqrt {s(x)}}}~\exp \left(-i\omega ~\int _{\zeta }^{x}s(y)~{\text{d}}y\right)~.$

If there are no sources in region III the solution decays exponentially in the $-x$ direction. Then the WKB solution with $s(x)=i\alpha (x)$ is

${\text{(34)}}\qquad \varphi _{III}(x)\sim {\cfrac {B_{-}}{\sqrt {\alpha (x)}}}~\exp \left(\omega ~\int _{\zeta }^{x}\alpha (y)~{\text{d}}y\right)$

where the coefficient $B_{-}=A_{-}/{\sqrt {i}}$.

In region II, since $\epsilon$ or $\mu$ vary linearly with $x$, we may write

${\text{(35)}}\qquad k_{3}^{2}=\omega ^{2}~\epsilon ~\mu -k_{2}^{2}\sim \Omega (x-\zeta )$

Then, from (7)

$\varphi _{II}(x)\sim C~{\text{Ai}}(-\eta )\qquad {\text{with}}\qquad \eta =\Omega ^{1/3}(x-\zeta )~.$

When $\omega$ is high, the region I, II, and III overlap. Also, from (35) we observe that $\Omega \propto \omega ^{2}$. Hence, the large $\eta$ expansion (equation (8)) for the Airy function can be used in the overlap region, i.e.,

$\varphi _{II}(x)\sim C~\pi ^{-1/2}~\eta ^{-1/4}~\sin \left({\cfrac {2}{3}}~\eta ^{3/2}+{\cfrac {\pi }{4}}\right)~.$

Substituting for $\eta$ and using the identity

$\sin \theta ={\cfrac {1}{2i}}\left(e^{i\theta }-e^{-i\theta }\right)$

we get

${\text{(36)}}\qquad \varphi _{II}(x)\sim {\cfrac {C}{2i~\pi ^{1/2}~\Omega ^{1/12}~(x-\zeta )^{1/4}}}~\left\{\exp \left[{\cfrac {2i}{3}}~\Omega ^{1/2}~(x-\zeta )^{3/2}+{\cfrac {\pi ~i}{4}}\right]-\exp \left[-{\cfrac {2i}{3}}~\Omega ^{1/2}~(x-\zeta )^{3/2}-{\cfrac {\pi ~i}{4}}\right]\right\}~.$

Also, in the neighborhood of region II,

$\omega ~s(x)\sim \Omega ^{1/2}~(x-\zeta )^{1/2}~.$

So

$\omega ~\int _{\zeta }^{x}s(y)~{\text{d}}y\sim {\cfrac {2}{3}}~\Omega ^{1/2}~(x-\zeta )^{3/2}~,~~x\sim \zeta ~.$

Therefore, $\varphi _{I}$ becomes

${\text{(37)}}\qquad \varphi _{I}(x)={\cfrac {A_{+}~\omega ^{1/2}}{\Omega ^{1/4}~(x-\zeta )^{1/4}}}~\exp \left({\cfrac {2i}{3}}~\Omega ^{1/2}~(x-\zeta )^{3/2}\right)+{\cfrac {A_{-}~\omega ^{1/2}}{\Omega ^{1/4}~(x-\zeta )^{1/4}}}~\exp \left(-{\cfrac {2i}{3}}~\Omega ^{1/2}~(x-\zeta )^{3/2}\right)~.$

Comparing (37) with (36) we get

${\cfrac {A_{+}}{A_{-}}}=-i\qquad {\text{and}}\qquad C=2~A_{+}~\omega ^{1/2}~\pi ^{1/2}~\Omega ^{-1/6}~e^{i\pi /4}~.$

Similarly, by comparing $\varphi _{II}$ and $\varphi _{III}$ in the region of overlap, we get

$C=2~B_{-}~\omega ^{1/2}~\pi ^{1/2}~\Omega ^{-1/6}~.$

References

• M. Abramowitz and I. A. Stegun. Airy functions. In Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, pages 446--452. Dover, New York, 1972.
• W. C. Chew. Waves and fields in inhomogeneous media. IEEE Press, New York, 1995.