Waves in composites and metamaterials/Airy solution and WKB solution

From Wikiversity

Jump to: navigation, search

Contents

[edit] 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


  \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 }{d x_3}\left(\cfrac{1}{\mu(x_3)}~\cfrac{d }{d x_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 k3 may be complex. Also, at the boundary, both \tilde{E}_1 and 1/\mu \partial\tilde{E}_1/\partial x_3 must be continuous.

[edit] TE waves in a non-magnetic medium

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

 \text{(4)} \qquad 
  \left[
  \cfrac{d^2 }{d x^2}
  + \omega^2~\epsilon~\mu - k_2^2\right]~\tilde{E}_1 = 0 \qquad 
  \text{where} \quad x \equiv x_3 ~.

[edit] 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 }{d x^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 Ai and 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 − η is shown in Figure~1.

Figure 1. The Airy functions Ai and Bi as functions of − η.

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

\begin{align}
  \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{align}

So Ai( − η) corresponds to an exponentially decaying wave as |x| \rightarrow \infty and Bi( − η) 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 Bi( − η) is unphysical which implies that C2 = 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 Ai( − η) 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.)

[edit] The Wentzel-Kramers-Brillouin (WKB) method

If we don't assume any particular linear variation of the permittivity ε(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}{d x^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}{d x^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 }{d x}\left(\epsilon^{-1}~\cfrac{d \tilde{H}_1}{d x}\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 }{d x}\left(\epsilon^{-1}~\cfrac{d \psi}{d x}\right) + 
  k_3^2(x)~\psi = 0 ~.

After expanding (11) we get

 \text{(12)} \qquad 
  \cfrac{d^2 \psi}{d x^2} = \epsilon^{-1}~\cfrac{d \epsilon}{d x}~\cfrac{d \psi}{d x}
  - 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}{d x^2} = \cfrac{d^2 }{d x^2}\left(\epsilon^{-1/2}\right)~\psi +
    \epsilon^{-1/2}~\cfrac{d^2 \psi}{d x^2} + 
    2~\cfrac{d }{d x}\left(\epsilon^{-1/2}\right)~\cfrac{d \psi}{d x} ~.

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

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

or,

 \text{(16)} \qquad 
  \cfrac{d^2 \varphi}{d x^2} + \left[k_3^2 - \epsilon^{1/2}~
    \cfrac{d^2 }{d x^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 k2 is proportional to ω which implies that k3 is also proportional to ω, i.e.,

 \text{(17)} \qquad 
  k_3^2(x) = \omega^2~s^2(x)

where s(x) is independent of ω.

In equation (16), if ω is large, then k3 will dominate and we will end up with exactly the same equation as (9), provided variations in ε are smooth (and we don't get large jumps in its derivatives).

Let us now try to solve (9). When k3 is constant, the solution of the equation is a traveling wave. If we assume that k3 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 τ(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., k3 is real) we can simplify the analysis slightly at this stage (even though this is not strictly necessary).

For large ω, 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 ω, 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 x0 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 ω 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 ω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 ω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 }{d x}\left[\ln(\omega~s(x))\right]

or,


  k_3(x) \gg \cfrac{d }{d x}\left[\ln(k_3(x))\right] ~.

Therefore, the restriction is that ω is large and that k3 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 k3(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 ω.

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 ζ 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α(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 ε or μ 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 ω is high, the region I, II, and III overlap. Also, from (35) we observe that \Omega \propto \omega^2. Hence, the large η 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 η 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}~.

[edit] Footnotes

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

[edit] 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.