# 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. [1] Then we may write

${\displaystyle {\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

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

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

${\displaystyle {\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

${\displaystyle 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:

${\displaystyle {\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

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

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

## TE waves in a non-magnetic medium

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

${\displaystyle {\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 ${\displaystyle x}$, then we may write

${\displaystyle \epsilon (x)=a+b~x}$

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

${\displaystyle {\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 ${\displaystyle B>0}$ (this is not strictly necessary, but simplifies things for our present analysis). Let us introduce a change of variables

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

Then (5) becomes

${\displaystyle {\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

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

where ${\displaystyle {\text{Ai}}}$ and ${\displaystyle {\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 ${\displaystyle -\eta }$ is shown in Figure~1.

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

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

{\displaystyle {\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 ${\displaystyle {\text{Ai}}(-\eta )}$ corresponds to an exponentially decaying wave as ${\displaystyle |x|\rightarrow \infty }$ and ${\displaystyle {\text{Bi}}(-\eta )}$ corresponds to an exponentially increasing waves at ${\displaystyle |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 ${\displaystyle x<0}$ then the solution ${\displaystyle {\text{Bi}}(-\eta )}$ is unphysical which implies that ${\displaystyle C_{2}=0}$. Therefore,

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

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

${\displaystyle {\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 ${\displaystyle \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

${\displaystyle {\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

${\displaystyle {\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

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

Recall also that the TM equation is

${\displaystyle {\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 ${\displaystyle \psi ={\tilde {H}}_{1}}$ to get

${\displaystyle {\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

${\displaystyle {\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

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

Differentiating (13) twice, we get

${\displaystyle {\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

${\displaystyle {\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,

${\displaystyle {\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

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

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

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

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

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

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

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

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

${\displaystyle {\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

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

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

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

${\displaystyle {\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

${\displaystyle {\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,

${\displaystyle {\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 ${\displaystyle \omega }$, equation (23) reduces to

${\displaystyle {\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,

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

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

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

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

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

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

Substituting (25) into (27) we get

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

or,

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

Integrating (28) we get

${\displaystyle {\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 ${\displaystyle \omega ^{2}}$ and higher) we get

${\displaystyle {\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

${\displaystyle {\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 ${\displaystyle k_{3}^{2}>0}$. Note that when ${\displaystyle s(x)=0}$, a solution does not exist.

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

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

Therefore,

${\displaystyle \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,

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

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

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

 Figure 3. Regions of validity of the linear solution and the WKB solution for large ${\displaystyle \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 ${\displaystyle \zeta }$ be the point on the ${\displaystyle x}$-axis where ${\displaystyle s(x)=0}$. Then, in region I, the solution is

${\displaystyle {\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 ${\displaystyle -x}$ direction. Then the WKB solution with ${\displaystyle s(x)=i\alpha (x)}$ is

${\displaystyle {\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 ${\displaystyle B_{-}=A_{-}/{\sqrt {i}}}$.

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

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

Then, from (7)

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

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

${\displaystyle \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 ${\displaystyle \eta }$ and using the identity

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

we get

${\displaystyle {\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,

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

So

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

Therefore, ${\displaystyle \varphi _{I}}$ becomes

${\displaystyle {\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

${\displaystyle {\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 ${\displaystyle \varphi _{II}}$ and ${\displaystyle \varphi _{III}}$ in the region of overlap, we get

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

## Footnotes

1. This content is based on the exposition in Chew95. Please consult that text for further details.

## 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.