Waves in composites and metamaterials/Transformation-based cloaking in electromagnetism

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
In this lecture we will give a brief description of cloaking in the context of conductivity. It is useful to start off with a desciption of some variational principles for electrical conductivity at this stage.

Variational principle
Suppose that the electrical conductivity $$\boldsymbol{\sigma}(\mathbf{x})$$ is real and symmetric. Also assume that

\alpha~\boldsymbol{\mathit{1}} \ge \boldsymbol{\sigma}(\mathbf{x}) \ge \beta~\boldsymbol{\mathit{1}} \qquad\text{for all}\alpha,\beta > 0 ~. $$ Consider the body ($$\Omega$$) with boundary ($$\partial\Omega$$) shown in Figure 1. We would like to minimize the power dissipation into heat inside the body. This statement can be expressed as

\min_{u ~,~ u=u_0~\rm{on}~\partial\Omega} W(u) $$ where

W(u) = \int_{\Omega} \boldsymbol{\nabla} u\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} u~\text{d}\Omega ~. $$ Now consider a variation $$v$$ where $$v = 0$$ on $$\partial\Omega$$ and let $$\delta$$ be a small parameter. Then

\begin{align} W(u + \delta~v) & = \int_{\Omega} \boldsymbol{\nabla} (u+\delta~v)\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} (u+\delta~v)~\text{d}\Omega \\ & =   \int_{\Omega} \boldsymbol{\nabla} u\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} u~\text{d}\Omega + 2~\delta~\int_{\Omega} \boldsymbol{\nabla} v\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} u~\text{d}\Omega + \delta^2~\int_{\Omega} \boldsymbol{\nabla} v\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} v~\text{d}\Omega > 0 ~. \end{align} $$ Using the identity

\mathbf{a}\cdot\boldsymbol{\nabla} b = \boldsymbol{\nabla} \cdot (b~\mathbf{a}) - b~\boldsymbol{\nabla} \cdot \mathbf{a} $$ in the middle term on the right hand side leads to

W(u + \delta~v) = \int_{\Omega} \boldsymbol{\nabla} u\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} u~\text{d}\Omega + 2~\delta~\int_{\Omega} \boldsymbol{\nabla} \cdot (v~\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)}~\text{d}\Omega - 2~\delta~\int_{\Omega} v~\boldsymbol{\nabla} \cdot (\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)}~\text{d}\Omega + \delta^2~\int_{\Omega} \boldsymbol{\nabla} v\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} v~\text{d}\Omega ~. $$ From the divergence theorem, we have

\int_{\Omega} \boldsymbol{\nabla} \cdot (v~\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)}~\text{d}\Omega = \int_{\partial\Omega} \mathbf{n}\cdot(v~\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} u)~\text{d}\Gamma $$ where $$\mathbf{n}$$ is the outward unit normal to the surface $$\partial\Omega$$ and $$\Gamma \equiv \partial\Omega$$. Since $$v = 0$$ on $$\partial\Omega$$, we have

\int_{\Omega} \boldsymbol{\nabla} \cdot (v~\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)}~\text{d}\Omega = 0 ~. $$ Therefore,

W(u + \delta~v) = \int_{\Omega} \boldsymbol{\nabla} u\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} u~\text{d}\Omega - 2~\delta~\int_{\Omega} v~\boldsymbol{\nabla} \cdot (\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)}~\text{d}\Omega + \delta^2~\int_{\Omega} \boldsymbol{\nabla} v\cdot\boldsymbol{\sigma}\cdot\boldsymbol{\nabla} v~\text{d}\Omega ~. $$ For $$W(u+\delta v)$$ to be positive for all $$v$$, it is sufficient to have

\int_{\Omega} v~\boldsymbol{\nabla} \cdot (\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)}~\text{d}\Omega = 0 ~. $$ If this is to be true for all $$v$$, then

{  \boldsymbol{\nabla} \cdot (\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)} = 0 ~. } $$ If we define the flux as

\mathbf{J}(\mathbf{x}) := \boldsymbol{\sigma}(\mathbf{x})\cdot\boldsymbol{\nabla} u $$ then we have

{ \boldsymbol{\nabla} \cdot \mathbf{J}(\mathbf{x}) = 0 ~. } $$

Coordinate transformation equations for currents
Let us take new curvilinear coordinates $$\mathbf{x}'(\mathbf{x})$$ as shown in Figure 2. The new coordinates are material coordinates.

The Jacobian of the transformation $$\mathbf{x}' \rightarrow \mathbf{x}$$ is given by

J = \det(\boldsymbol{A}) ~; A_{ij} := \frac{\partial x'_i}{\partial x_j} ~. $$ Then an infinitesimal volume $$\text{d}\Omega$$ of the body transforms as

\text{d}\Omega' \rightarrow J~\text{d}\Omega ~. $$ Recall that

W(u) = \int_{\Omega} \frac{\partial u}{\partial x_i}~\sigma_{ij}~\frac{\partial u}{\partial x_j}~\text{d}\Omega ~. $$ Then, using the chain rule, we get

W(u) = \int_{\Omega'} \left(\frac{\partial u}{\partial x'_m}~\frac{\partial x'_m}{\partial x_i}\right) ~\sigma_{ij}~\left(\frac{\partial x'_l}{\partial x_j}~\frac{\partial u}{\partial x'_l}\right) ~\cfrac{1}{J}~\text{d}\Omega' $$ or,

W(u) = \int_{\Omega'} \frac{\partial u}{\partial x'_m}~\sigma'_{ml}~\frac{\partial u}{\partial x'_l} ~\text{d}\Omega' $$ where

\sigma'_{ml} = \cfrac{1}{J}~\frac{\partial x'_m}{\partial x_i}~\sigma_{ij}~ \frac{\partial x'_l}{\partial x_j} = \cfrac{1}{J}~A_{mi}~\sigma_{ij}~A_{lj} ~. $$ Hence, in the transformed coordinates, the functional $$W(u)$$ takes the form

W(u) = \int_{\Omega'} \boldsymbol{\nabla}' u\cdot\boldsymbol{\sigma}'\cdot\boldsymbol{\nabla}' u~\text{d}\Omega' $$ where $$\boldsymbol{\nabla}' (\bullet)$$ denotes a gradient with respect to the $$\mathbf{x}'$$ coordinates and the conductivity transforms as

\boldsymbol{\sigma}'(\mathbf{x}') = \cfrac{1}{J}~\boldsymbol{A}(\mathbf{x})\cdot\boldsymbol{\sigma}(\mathbf{x})\cdot\boldsymbol{A}^T(\mathbf{x}) ~. $$

Interpretation
We can now interpret the minimization problem in the transformed coordinates as follows:
 * The function $$u'(x') = u(x')$$ minimizes $$W$$ in a body $$\Omega'$$     filled with material with conductivity $$\boldsymbol{\sigma}'(\mathbf{x}')$$ with $$x'_1,x'_2,x'_3$$      as Cartesian coordinates in $$x'$$ space.

Therefore, for $$W$$ to remain positive, we must have

\mathbf{J}'(\mathbf{x}') = \boldsymbol{\sigma}'(\mathbf{x}') \cdot \boldsymbol{\nabla}' u(\mathbf{x}') = \cfrac{1}{J}~[\boldsymbol{A}(\mathbf{x})\cdot\boldsymbol{\sigma}(\mathbf{x})\cdot\boldsymbol{A}^T(\mathbf{x})] \cdot \boldsymbol{\nabla}' u(\mathbf{x}')~. $$ Now,

[\boldsymbol{\nabla} u]_i = \frac{\partial u}{\partial x_i} = \frac{\partial x'_m}{\partial x_i}~\frac{\partial u}{\partial x'_m} = A_{mi}~\frac{\partial u}{\partial x'_m} = [\boldsymbol{A}^T\cdot\boldsymbol{\nabla}' u]_i ~. $$ Hence,

\mathbf{J}'(\mathbf{x}') = \cfrac{1}{J}~[\boldsymbol{A}(\mathbf{x})\cdot\boldsymbol{\sigma}(\mathbf{x})\cdot\boldsymbol{A}^T(\mathbf{x})] \cdot [\boldsymbol{A}^T(\mathbf{x})]^{-1}\cdot\boldsymbol{\nabla} u(\mathbf{x}) $$ or,

{ \mathbf{J}'(\mathbf{x}') = \cfrac{1}{J}~\boldsymbol{A}(\mathbf{x})\cdot\boldsymbol{\sigma}(\mathbf{x})\cdot\boldsymbol{\nabla} u(\mathbf{x}) = \cfrac{\boldsymbol{A}\cdot\mathbf{J}(\mathbf{x})}{\det(\boldsymbol{A})} ~. } $$ This is the transformation law for currents. Using the same arguments as before, we can show that

{ \boldsymbol{\nabla}'\cdot\mathbf{J}'(\mathbf{x}') = 0 ~. } $$ Let the electric field $$\mathbf{E}$$ be derived from the potential $$u$$. Then the fields

\mathbf{E}(\mathbf{x}) = \boldsymbol{\nabla} u(\mathbf{x}) \qquad \text{and} \qquad \mathbf{E}'(\mathbf{x}') = \boldsymbol{\nabla}' u(\mathbf{x}') $$ are related via

{  \mathbf{E}'(\mathbf{x}') = (\boldsymbol{A}^T)^{-1}\cdot\mathbf{E}(\mathbf{x}) ~. } $$ Therefore, there are two transformations which are equivalent. However, an isotropic material transforms to an anisotropic material via the transformation equation for conductivity.

Electrical tomography
Consider the situation shown in Figure 1. Let the conductivity of the body be $$\boldsymbol{\sigma}(\mathbf{x})$$ and let us require that $$\boldsymbol{\nabla} \cdot (\boldsymbol{\sigma}\cdot\boldsymbol{\nabla u)} = 0$$ inside the body. In electrical tomography one measures the current flux $$\mathbf{n}\cdot\mathbf{J}(\mathbf{x})$$ at the surface for all choices of the potential $$u_0$$.

Suppose one knows the Dirchlet to Neumann map ($$\varphi$$)

\varphi : u_0 \rightarrow \mathbf{n}\cdot\mathbf{J}(\mathbf{x}) = g(\mathbf{x}) ~. $$ Can one find $$\boldsymbol{\sigma}(\mathbf{x})$$? No, not generally. Figure 3 illustrates why that is the case. For the body in the figure, the transformation is $$\mathbf{x}' = \mathbf{x}$$ outside the blue region while inside the blue region $$\mathbf{x}' \ne \mathbf{x}$$. Also, outside the blue region, $$\Omega' = \Omega$$, $$\mathbf{J}' = \mathbf{J}$$, and $$u' = u$$. Inside the blue region $$\Omega' \ne \Omega$$ and $$\mathbf{J}'$$ is obtained via the transformation rule.

From the figure we can see that the Dirichlet-Neumman map will remain unchanged on $$\partial\Omega$$. Hence, the body appears to be exactly the same in $$\mathbf{x}'$$-space but has a different conductivity.

Even though this fact has been known for a while, there was still hope that you could determine $$\boldsymbol{\sigma}(\mathbf{x})$$ uniquely, modulo a coordinate transformation. However, such hopes were dashed when Greenleaf, Lassas, and Uhlmann provided a counterexample in 2003 (Greenleaf03).

First transformation based example of cloaking
Greenleaf et al. (Greenleaf03) provided the first example of transformation based cloaking. They considered a singular transformation

\mathbf{x}'(\mathbf{x}) = \begin{cases} \left(\cfrac{|\mathbf{x}|}{2} + 1\right)~\cfrac{\mathbf{x}}{|\mathbf{x}|} & \text{if}~|\mathbf{x}| < 2 \\ \mathbf{x} & \text{if}~|\mathbf{x}| > 2 ~. \end{cases} $$ The effect of this mapping is shown in the schematic in Figure 4. An epsilon ball at the center of $$\Omega$$ is mapped into a sphere of radius 1 in $$\Omega'$$. The value of $$\boldsymbol{\sigma}'(\mathbf{x}')$$ is singular at the boundary of this sphere. Inside the sphere of radius 1, the transformed conductivity has the form $$\boldsymbol{\sigma}'(\mathbf{x}') = h(\mathbf{x}')$$.

Therefore we can put a small body inside and the potential outside will be undisturbed by the presence of the body in the cloaking region.

Cloaking for Electromagnetism
Pendry, Schurig, and Smith (Pendry06) showed in 2006 that cloaking could be achieved for electromagnetic waves. The concept of cloaking follows from the observation that Maxwell's equations keep their form under coordinate transformations. The Maxwell's equations at fixed frequency $$\omega$$ are

\begin{align} \boldsymbol{\nabla} \times \mathbf{E} + i\omega\boldsymbol{\mu}\cdot\mathbf{H} &= \boldsymbol{0} \\ \boldsymbol{\nabla} \times \mathbf{H} - i\omega\boldsymbol{\epsilon}\cdot\mathbf{E} &= \boldsymbol{0} ~. \end{align} $$ A coordinate transformation ($$\mathbf{x} \rightarrow \mathbf{x}'$$) gives us the equivalent relations

\begin{align} \boldsymbol{\nabla}'\times\mathbf{E}' + i\omega\boldsymbol{\mu}'\cdot\mathbf{H}' &= \boldsymbol{0} \\ \boldsymbol{\nabla}'\times\mathbf{H}' - i\omega\boldsymbol{\epsilon}'\cdot\mathbf{E}' &= \boldsymbol{0} \end{align} $$ with

\mathbf{x}' = \mathbf{x}'(\mathbf{x}) ~; \mathbf{E}'(\mathbf{x}') = (\boldsymbol{A}^T)^{-1}~\mathbf{E}(\mathbf{x}) ~; \mathbf{H}'(\mathbf{x}') = (\boldsymbol{A}^T)^{-1}~\mathbf{H}(\mathbf{x}) $$ where

A_{ki} = \frac{\partial x'_k}{\partial x_i} ~; [\boldsymbol{A}^{-1}]_{ij} = \frac{\partial x_i}{\partial x'_j} $$ and

\boldsymbol{\mu}'(\mathbf{x}') = \cfrac{\boldsymbol{A}\cdot\boldsymbol{\mu}(\mathbf{x})\cdot\boldsymbol{A}^T}{\det(\boldsymbol{A})} ~; \boldsymbol{\epsilon}'(\mathbf{x}') = \cfrac{\boldsymbol{A}\cdot\boldsymbol{\epsilon}(\mathbf{x})\cdot\boldsymbol{A}^T}{\det(\boldsymbol{A})} ~. $$ To see that this invariance of form under coordinate transformations does indeed hold, observe that
 * $$ \text{(1)} \qquad

-i\omega\boldsymbol{\mu}'\cdot\mathbf{H}' = \cfrac{-i\omega\boldsymbol{A}\cdot\boldsymbol{\mu}\cdot\mathbf{H}}{\det(\boldsymbol{A})} = \cfrac{\boldsymbol{A}\cdot(\boldsymbol{\nabla} \times \mathbf{E})}{\det(\boldsymbol{A})} ~. $$ We want to show that this equals $$\boldsymbol{\nabla}'\times\mathbf{E}'$$.

In index notation, (1) can be written as

\begin{align} \left[-i\omega\boldsymbol{\mu}^{'}\cdot\mathbf{H}^{'}\right]_h & = \cfrac{\left[\boldsymbol{A}\cdot\left(\boldsymbol{\nabla} \times \mathbf{E}\right)\right]_h}{\det(\boldsymbol{A})} \\ & = \cfrac{1}{\det(\boldsymbol{A})}~\frac{\partial x'_h}{\partial x_j}~\mathcal{E}_{jmk}~\frac{\partial E_k}{\partial x_m} \\ & = \cfrac{1}{\det(\boldsymbol{A})}~\frac{\partial x'_h}{\partial x_j}~\mathcal{E}_{jmk}~\frac{\partial x'_l}{\partial x_m} ~\frac{\partial E_k}{\partial x'_l} ~. \end{align} $$ On the other hand,

\begin{align} \left[\boldsymbol{\nabla}'\times\mathbf{E}'\right]_h & = \left[\boldsymbol{\nabla}'\times(\boldsymbol{A}^{-T}\cdot\mathbf{E})\right]_h \\ & = \mathcal{E}_{hlm}~\frac{\partial }{\partial x_i}\left(\frac{\partial x_k}{\partial x'_m}~E_k\right) \\ & = \mathcal{E}_{hlm}~\frac{\partial^2 x_k}{\partial x'_l \partial x'_m}~E_k + \mathcal{E}_{hlm}~\frac{\partial x_k}{\partial x'_m}~\frac{\partial E_k}{\partial x'_l} \\ & = \mathcal{E}_{hlm}~\frac{\partial x_k}{\partial x'_m}~\frac{\partial E_k}{\partial x'_l} ~. \end{align} $$ The first term above evaluates to zero because of $$\text{tr}(\boldsymbol{A}\cdot\boldsymbol{B}) = 0$$ if $$\boldsymbol{A}$$ is skew and $$\boldsymbol{B}$$ is symmetric.

So we now need to show that

\mathcal{E}_{hlm}~\frac{\partial x_k}{\partial x'_m}~\frac{\partial E_k}{\partial x'_l} = \cfrac{1}{\det(\boldsymbol{A})}~\frac{\partial x'_h}{\partial x_j}~\mathcal{E}_{jmk}~\frac{\partial x'_l}{\partial x_m} ~\frac{\partial E_k}{\partial x'_l} $$ or that,
 * $$ \text{(2)} \qquad

\det(\boldsymbol{A})~\mathcal{E}_{hlm}~\frac{\partial x_k}{\partial x'_m} = \frac{\partial x'_h}{\partial x_j}~\mathcal{E}_{jmk}~\frac{\partial x'_l}{\partial x_m} ~. $$ Multiply both sides of (2) by $$A_{pk}$$ and sum over $$k$$, (i.e., multiply by $$\boldsymbol{A}$$ which is non-singular). Then we get

\det(\boldsymbol{A})~\mathcal{E}_{hlm}~\frac{\partial x_k}{\partial x'_m}~\frac{\partial x'_p}{\partial x_k} = \mathcal{E}_{jmk}~\frac{\partial x'_h}{\partial x_j}~\frac{\partial x'_l}{\partial x_m}~\frac{\partial x'_p}{\partial x_k} $$ or,

\det(\boldsymbol{A})~\mathcal{E}_{hlm}~\frac{\partial x'_p}{\partial x'_m} = \det(\boldsymbol{A})~\mathcal{E}_{hlm}~\delta_{pm} = \det(\boldsymbol{A})~\mathcal{E}_{hlp} = \mathcal{E}_{jmk}~\frac{\partial x'_h}{\partial x_j}~\frac{\partial x'_l}{\partial x_m}~\frac{\partial x'_p}{\partial x_k} ~. $$ Both sides are completely antisymmetric with respect o $$h,l,p$$. So it suffices to take $$h = 1$$, $$l = 2$$, $$p = 3$$ and we can write

\det(\boldsymbol{A})~\mathcal{E}_{123} = \det(\boldsymbol{A}) = \mathcal{E}_{jmk}~\frac{\partial x'_1}{\partial x_j}~\frac{\partial x'_2}{\partial x_m}~\frac{\partial x'_3}{\partial x_k} ~. $$ The right hand side above is the well known formula for the determinant of the Jacobian. Hence the first of the transformed Maxwell equations holds. We can follow the same procedure to show that the second Maxwell's equation also maintains its form under coordinate transformations. Hence Maxwell's equations are invariant with respect to coordinate transformations.