Waves in composites and metamaterials/Effective tensors using Hilbert space formalism

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.

Recap

In the previous lecture we introduced the Hilbert space $\mathscr{H}$ of of square-integrable, $Q$-periodic, complex vector fields with the inner product

$(𝐚,𝐛)=\underset{Q}{\int }\overline{𝐚(𝐱)}\cdot 𝐛(𝐱)\text{d}𝐱\equiv \underset{\hat{Q}}{\int }\widehat{𝐚}(𝐤)\cdot \widehat{𝐛}(𝐤)\text{d}𝐤.$

We then decomposed $\mathscr{H}$ into three orthogonal subspaces.

1. The subspace $𝒰$ of uniform fields, i.e., $𝐚(𝐱)$ is independent of $𝐱$, or in Fourier space, $\widehat{𝐚}(𝐤)=0$ unless $𝐤=0$.
2. The subspace $𝒥$ of zero divergence, zero average value fields, i.e., $\mathrm{\nabla }\cdot 𝐚(𝐱)=0$ and $⟨𝐚(𝐱)⟩=0$, or in Fourier space, $\widehat{𝐚}(0)=0$ and $𝐤\cdot \widehat{𝐚}(𝐤)=0$.
3. The subspace $ℰ$ of zero curl, zero average value fields, i.e., $\mathrm{\nabla }\times 𝐚(𝐱)=0$ and $⟨𝐚(𝐱)⟩=0$, or in Fourier space, $\widehat{𝐚}(0)=0$ and $\widehat{𝐚}(𝐤)=\alpha (𝐤)𝐤$.

To determine the effective permittivity, we introduced the operator $\mathrm{\Gamma }$ with the properties

$\text{(1)}𝐛(𝐱)=\mathrm{\Gamma }\cdot 𝐚(𝐱)$

if and only if

$𝐛\in ℰ\text{and}𝐚-{\mathit{ϵ}}_0\cdot 𝐛\in 𝒰\oplus 𝒥⟹\mathrm{\nabla }\cdot (𝐚-{\mathit{ϵ}}_0\cdot 𝐛)=0.$

We also found that in Fourier space,

$\widehat{𝐛}(𝐤)=\hat{\mathrm{\Gamma }}(𝐤)\cdot \widehat{𝐚}(𝐤)$

where

$\text{(2)}\hat{\mathrm{\Gamma }}(𝐤)=\{\begin{array}{cc}\frac{𝐤\otimes 𝐤}{𝐤\cdot {\mathit{ϵ}}_0\cdot 𝐤}=\mathrm{\Gamma }(𝐧)& \text{with}𝐧=\frac{𝐤}{|𝐤|}\text{if}𝐤\ne 0\\ 0& \text{if}𝐤=0.\end{array}$

Deriving a formula for the effective permittivity

Let us now derive a formula for the effective tensor. Recall that the polarization is defined as

$\text{(3)}𝐏=(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot 𝐄=𝐃-{\mathit{ϵ}}_0\cdot 𝐄$

where the permittivity tensor is being thought of as an operator that operates on $𝐄$.

Also notice that

$\text{(4)}⟨𝐄⟩-𝐄\in ℰ\text{(curl free)}.$

From (3) we have

$\text{(5)}𝐏-{\mathit{ϵ}}_0\cdot [⟨𝐄⟩-𝐄]=𝐃-{\mathit{ϵ}}_0\cdot ⟨𝐄⟩\in 𝒰\oplus 𝒥\text{(divergence free)}.$

From the definition of $\mathrm{\Gamma }$ (equations (1) and (2)) and using (4) and (5), we can show that

$\text{(6)}\mathrm{\Gamma }\cdot 𝐏=⟨𝐄⟩-𝐄.$

Let $\mathit{ϵ}-{\mathit{ϵ}}_0$ act on both sides of (6). Then we get

$(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot \mathrm{\Gamma }\cdot 𝐏=(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot ⟨𝐄⟩-(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot 𝐄=(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot ⟨𝐄⟩-𝐏.$

Therefore,

$\text{(7)}[1+(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot \mathrm{\Gamma }]\cdot 𝐏=(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot ⟨𝐄⟩.$

Inverting (7) gives

$\text{(8)}𝐏=[1+(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot \mathrm{\Gamma }{]}^{-1}\cdot (\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot ⟨𝐄⟩.$

Averaging (8) gives

$\text{(9)}⟨𝐏⟩=⟨[1+(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot \mathrm{\Gamma }{]}^{-1}\cdot (\mathit{ϵ}-{\mathit{ϵ}}_0)⟩\cdot ⟨𝐄⟩.$

Averaging (3) leads to the relation

$\text{(10)}⟨𝐏⟩=({\mathit{ϵ}}_{\text{eff}}-{\mathit{ϵ}}_0)\cdot ⟨𝐄⟩.$

Comparing (9) and (10) shows us that

$\text{(11)}{\mathit{ϵ}}_{\text{eff}}={\mathit{ϵ}}_0+⟨[1+(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot \mathrm{\Gamma }{]}^{-1}\cdot (\mathit{ϵ}-{\mathit{ϵ}}_0)⟩.$

Recall from the previous lecture that, in equation (11), the operator $\mathit{ϵ}$ is local in real space while the operator $\mathrm{\Gamma }$ is local in Fourier space. It is therefore not obvious how one can invert $[1+(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot \mathrm{\Gamma }]$.

Let us define

$\mathit{\delta }:=\mathit{ϵ}-{\mathit{ϵ}}_0.$

Then (8) can be written as

$𝐏=[1+\mathit{\delta }\cdot \mathrm{\Gamma }{]}^{-1}\cdot \mathit{\delta }\cdot ⟨𝐄⟩.$

Assuming that $-\mathit{\delta }\cdot \mathrm{\Gamma }<1$, we can expand the first operator in terms of an infinite series, i.e.,

$[1+\mathit{\delta }\cdot \mathrm{\Gamma }{]}^{-1}=\sum _{n=0}^{\mathrm{\infty }}(-\mathit{\delta }\cdot \mathrm{\Gamma }{)}^n.$

Then we have

$𝐏=[\sum _{n=0}^{\mathrm{\infty }}(-\mathit{\delta }\cdot \mathrm{\Gamma }{)}^n]\cdot \mathit{\delta }\cdot ⟨𝐄⟩.$

Also, from the definition of $𝐏$, we have

$𝐏=\mathit{\delta }\cdot 𝐄.$

Hence,

$\mathit{\delta }\cdot 𝐄=[\sum _{n=0}^{\mathrm{\infty }}(-\mathit{\delta }\cdot \mathrm{\Gamma }{)}^n]\cdot \mathit{\delta }\cdot ⟨𝐄⟩.$

Now,

$(-\mathit{\delta }\cdot \mathrm{\Gamma }{)}^n\cdot \mathit{\delta }=(-1{)}^n\mathit{\delta }\cdot \mathrm{\Gamma }\cdot \mathit{\delta }\cdot \mathrm{\Gamma }\dots \mathit{\delta }\cdot \mathrm{\Gamma }\cdot \mathit{\delta }=\mathit{\delta }\cdot (-\mathrm{\Gamma }\cdot \mathit{\delta }{)}^n.$

Therefore,

$\mathit{\delta }\cdot 𝐄=\mathit{\delta }\cdot [\sum _{n=0}^{\mathrm{\infty }}(-\mathrm{\Gamma }\cdot \mathit{\delta }{)}^n]\cdot ⟨𝐄⟩$

or,

$\text{(12)}𝐄=[\sum _{n=0}^{\mathrm{\infty }}(-\mathrm{\Gamma }\cdot \mathit{\delta }{)}^n]\cdot ⟨𝐄⟩.$

Let us now define

$\text{(13)}\begin{array}{cc}{𝐏}^{(m)}& :=[\sum _{n=0}^m(-\mathit{\delta }\cdot \mathrm{\Gamma }{)}^n]\cdot \mathit{\delta }\cdot ⟨𝐄⟩\\ {𝐄}^{(m)}& :=[\sum _{n=0}^m(-\mathrm{\Gamma }\cdot \mathit{\delta }{)}^n]\cdot ⟨𝐄⟩.\end{array}$

Then we can write

$\text{(14)}\mathrm{\Gamma }\cdot {𝐏}^{(m)}=⟨𝐄⟩-{𝐄}^{(m+1)}\text{and}\mathit{\delta }\cdot {𝐄}^{(m)}={𝐏}^{(m)}.$

These recurrence relations may be used to compute these fields inductively. An algorithm that may be used is outlined below:

• Set $m=0$. Then

${𝐄}^{(0)}=⟨𝐄⟩.$

• Compute ${𝐏}^{(m)}(𝐱)$ in real space using the relation

${𝐏}^{(m)}=\mathit{\delta }\cdot {𝐄}^{(m)}.$

• Take a fast Fourier transform to find ${\widehat{𝐏}}^{(m)}(𝐤)$.
• From (14) we get

${\widehat{𝐄}}^{(m+1)}=\{\begin{array}{cc}-\hat{\mathrm{\Gamma }}(𝐤)\cdot {\widehat{𝐏}}^{(m)}(𝐤)& \text{for}𝐤\ne 0\\ ⟨𝐄⟩& \text{for}𝐤=0.\end{array}$

• Compute ${\widehat{𝐄}}^{(m+1)}(𝐤)$ in Fourier space.
• Take an inverse fast Fourier transform to find ${𝐄}^{(m+1)}(𝐱)$ in real space.
• Increment $m$ by 1 and repeat.

This is the method of Moulinec and Suquet (Mouli94). The method also extends to nonlinear problems (Mouli98). However, there are other iterative methods that have faster convergence.

Convergence of expansions

For simplicity, let us assume that ${\mathit{ϵ}}_0$ is isotropic, i.e., ${\mathit{ϵ}}_0={ϵ}_{0}1$. Then,

$\mathrm{\Gamma }=\frac{{\mathrm{\Gamma }}_1}{{ϵ}_0}$

where ${\mathrm{\Gamma }}_1$ is the projection from $\mathscr{H}$ onto $ℰ$.

Define the norm of a field $𝐏$ as

$|𝐏|:=(𝐏,𝐏{)}^{1/2}=⟨\stackrel{1/2}{\overline{𝐏⟩\cdot 𝐏}}.$

Also, define the norm of a linear operator $𝑨$ acting on $𝐏$ as

$\Vert 𝑨{\Vert }_{}:=\underset{𝐏\in \mathscr{H},|𝐏|=1}{\sup }|𝑨\cdot 𝐏|.$

Therefore,

$|𝑨\cdot 𝐏|\le \Vert 𝑨{\Vert }_{}\cdot |𝐏|.$

Hence,

$|𝑨\cdot 𝑩\cdot 𝐏|\le \Vert 𝑨{\Vert }_{}\cdot |𝑩\cdot 𝐏|\le \Vert 𝑨{\Vert }_{}\cdot \Vert 𝑩\Vert \cdot |𝐏|.$

So

$\Vert 𝑨\cdot 𝑩{\Vert }_{}\le \Vert 𝑨{\Vert }_{}\cdot \Vert 𝑩{\Vert }_{}.$

In addition, we have the triangle inequality

$\Vert 𝑨+𝑩{\Vert }_{}\le \Vert 𝑨{\Vert }_{}+\Vert 𝑩{\Vert }_{}.$

So from (12) and (13), we have

$\begin{array}{cc}|𝐄-{𝐄}^{(m)}|& =|\sum _{n=m+1}^{\mathrm{\infty }}(-\mathit{\delta }\cdot \mathrm{\Gamma }{)}^n\cdot ⟨𝐄⟩|\\ & \le \Vert \sum _{n=m+1}^{\mathrm{\infty }}(-\mathit{\delta }\cdot \mathrm{\Gamma }{)}^n{\Vert }_{}\cdot |⟨𝐄⟩|\\ & \le (\sum _{n=m+1}^{\mathrm{\infty }}\Vert \frac{\mathit{\delta }}{{ϵ}_0}{\Vert }^n\cdot \Vert {\mathrm{\Gamma }}_1\underset{}{\overset{n}{\Vert }})\cdot |⟨𝐄⟩|.\end{array}$

But $\Vert {\mathrm{\Gamma }}_1{\Vert }_{}=1$ since it is a projection onto $ℰ$. Hence,

$|𝐄-{𝐄}^{(m)}|\le (\sum _{n=m+1}^{\mathrm{\infty }}\Vert \frac{\mathit{\delta }}{{ϵ}_0}{\Vert }^n)\cdot |⟨𝐄⟩|.$

Therefore the series converges provided

$\Vert \frac{\mathit{\delta }}{{ϵ}_0}\Vert <1.$

In that case

$|𝐄-{𝐄}^{(m)}|\le \frac{{\gamma }^{m+1}}{1-\gamma }|⟨𝐄⟩|$

where

$\gamma :=\Vert \frac{\mathit{\delta }}{{ϵ}_0}{\Vert }_{}.$

To get a better understanding of the norm $\gamma$, let us consider a $n$-phase composite with isotropic phases, i.e.,

$\mathit{ϵ}=\sum _{i=1}^n{\chi }_i{ϵ}_{i}1$

where

${\chi }_i=\{\begin{array}{cc}1& \text{in phase}i\\ 0& \text{otherwise}.\end{array}$

In this case,

$\gamma =\Vert \frac{\mathit{\delta }}{{ϵ}_0}\Vert =\underset{𝐏\in \mathscr{H},|𝐏|=1}{\sup }\sum _{i=1}^n{⟨\overline{\left(\frac{{ϵ}_i-{ϵ}_0}{{ϵ}_0}{𝐏}_i\right)}\cdot \left(\frac{{ϵ}_i-{ϵ}_0}{{ϵ}_0}{𝐏}_i\right)⟩}^{1/2}$

where

${𝐏}_i={\chi }_i𝐏\text{and}|𝐏|={(\sum _i\overline{{𝐏}_i}\cdot 𝐏)}^{1/2}=1.$

Hence,

$\gamma =\underset{𝐏\in \mathscr{H},|𝐏|=1}{\sup }\sum _{i=1}^n\left|\frac{{ϵ}_i-{ϵ}_0}{{ϵ}_0}\right|⟨{\overline{𝐏}}_i\cdot {𝐏}_i{⟩}^{1/2}.$

Since, $⟨{\overline{𝐏}}_i\cdot {𝐏}_i{⟩}^{1/2}$ are weights, it makes sense to put the weights where ${ϵ}_i$ is maximum. Hence, we can write

$\gamma =\underset{i}{\max }\left|\frac{{ϵ}_i-{ϵ}_0}{{ϵ}_0}\right|=\underset{i}{\max }\frac{|{ϵ}_i-{ϵ}_0|}{|{ϵ}_0|}.$

For $\gamma$ to be less than 1, we therefore require that, for all $i$,

$|{ϵ}_i-{ϵ}_0|<|{ϵ}_0|.$

A geometrical representation of this situation is shown in Figure 1.

File:Lec21Fig1.jpg Figure 1. Allowable values of ${ϵ}_i$ for convergence of series expansion.

If the value of ${ϵ}_0$ is sufficiently large, then we get convergence if all the ${ϵ}_i$ s lie in one half of the complex plane (shown by the green line in the figure).

Similarly, we can expand

${\mathit{ϵ}}_{\text{eff}}={\mathit{ϵ}}_0+⟨[1+(\mathit{ϵ}-{\mathit{ϵ}}_0)\cdot \mathrm{\Gamma }{]}^{-1}\cdot (\mathit{ϵ}-{\mathit{ϵ}}_0)⟩$

in the form

$\text{(15)}{\mathit{ϵ}}_{\text{eff}}={\mathit{ϵ}}_0+\sum _{j=0}^{\mathrm{\infty }}{\mathrm{\Gamma }}_0\cdot \mathit{\delta }\cdot (-\mathrm{\Gamma }\cdot \mathit{\delta }{)}^j\cdot {\mathrm{\Gamma }}_0$

where ${\mathrm{\Gamma }}_0$ is a projection onto $𝒰$, i.e.,

${\mathrm{\Gamma }}_0\cdot 𝐏=⟨𝐏⟩.$

In this case, we find that the series converges provided

$|{ϵ}_i-{ϵ}_0|<|{ϵ}_0|\text{for all}i.$

Note that each term in (15) is an analytic function of ${ϵ}_1$ (in fact a polynomial). So, if we truncate the series, we have an analytic function of ${ϵ}_1$.

Since a sequence of analytic functions which is uniformly convergent in a domain converges to a function which is analytic in that domain (see, for example, Rudin76), we deduce that ${\mathit{ϵ}}_{\text{eff}}$ is an analytic function of ${ϵ}_1$ in the ${ϵ}_0$ disk (see Figure~1) with $r<|{ϵ}_0|$ provided $|{ϵ}_i-{ϵ}_0|<{ϵ}_0$ for $i=2,3,\dots ,n$.

Similarly, the effective tensor is an analytic function of ${ϵ}_2$, ${ϵ}_3$, etc.

Since ${\mathit{ϵ}}_{\text{eff}}$ is independent of ${ϵ}_0$, by taking the union of all such regions of analyticity, we conclude that ${\mathit{ϵ}}_{\text{eff}}$ is an analytic function of ${ϵ}_1,\dots ,{ϵ}_n$ provided all these ${ϵ}_i$ s lie inside a half-plane (see Figure~1). This means that there exists a $\theta$ such that

$\text{Re}(e^{i\theta }{ϵ}_j)>0\text{for all}j.$

Corollary:

A corollary of the above observations is the following. If each ${ϵ}_i(\omega )$ is an analytic function of $\omega$ for $\text{Im}(\omega )>0$ (which is what one expects with $\omega$ as the frequency) and $\text{Im}[{ϵ}_i(\omega )]>0$ for all $\omega$ with $\text{Re}\omega >0$, then ${\mathit{ϵ}}_{\text{eff}}(\omega )$ will be analytic in $\omega$.

Another interesting property:

Now, if $𝐃=\mathit{ϵ}\cdot 𝐄$, we have

$\lambda 𝐃=\lambda \mathit{ϵ}\cdot 𝐄\text{for all constants}\lambda .$

Therefore,

$⟨𝐃⟩={\mathit{ϵ}}_{\text{eff}}\cdot ⟨𝐄⟩⟹⟨\lambda 𝐃⟩=\lambda {\mathit{ϵ}}_{\text{eff}}\cdot ⟨𝐄⟩.$

This means that

$(\lambda \mathit{ϵ}{)}_{\text{eff}}=\lambda {\mathit{ϵ}}_{\text{eff}}.$

Therefore, ${\mathit{ϵ}}_{\text{eff}}({ϵ}_1,{ϵ}_2,\dots ,{ϵ}_n)$ is homogeneous of degree one and

${\mathit{ϵ}}_{\text{eff}}(\lambda {ϵ}_1,\lambda {ϵ}_2,\dots ,\lambda {ϵ}_n)=\lambda {\mathit{ϵ}}_{\text{eff}}({ϵ}_1,{ϵ}_2,\dots ,{ϵ}_n).$

For a two-phase composite, if we take

$\lambda =\frac{1}{{ϵ}_2}$

we get

${\mathit{ϵ}}_{\text{eff}}({ϵ}_1,{ϵ}_2)={ϵ}_2{\mathit{ϵ}}_{\text{eff}}({ϵ}_1/{ϵ}_2,1).$

Therefore, it suffices to study the analytic function ${\mathit{ϵ}}_{\text{eff}}({ϵ}_1,1)={\mathit{ϵ}}_{\text{eff}}({ϵ}_1)$. For further details see Milton02.

References

• Template:Wikicite     G. W. Milton. Theory of Composites. Cambridge University Press, New York, 2002.
• Template:Wikicite     H. Moulinec and P. Suquet. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes rendus de l'Académie des sciences II, 318(11):1417--1423, 1994.
• Template:Wikicite     H. Moulinec and P. Suquet. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Methods Appl. Mech. Engrg., 157:69--94, 1998.
• Template:Wikicite     W. Rudin. Principles of Mathematical Analysis. McGraw-Hill, New York, 1976.