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.
In the previous lecture we introduced the Hilbert space
H
{\displaystyle {\mathcal {H}}}
of
of square-integrable,
Q
{\displaystyle Q}
-periodic, complex vector fields with the
inner product
(
a
,
b
)
=
∫
Q
a
(
x
)
¯
⋅
b
(
x
)
d
x
≡
∫
Q
^
a
^
(
k
)
⋅
b
^
(
k
)
d
k
.
{\displaystyle (\mathbf {a} ,\mathbf {b} )=\int _{Q}{\overline {\mathbf {a} (\mathbf {x} )}}\cdot \mathbf {b} (\mathbf {x} )~{\text{d}}\mathbf {x} \equiv \int _{\hat {Q}}{\widehat {\mathbf {a} }}(\mathbf {k} )\cdot {\widehat {\mathbf {b} }}(\mathbf {k} )~{\text{d}}\mathbf {k} ~.}
We then decomposed
H
{\displaystyle {\mathcal {H}}}
into three orthogonal subspaces.
The subspace
U
{\displaystyle {\mathcal {U}}}
of uniform fields, i.e.,
a
(
x
)
{\displaystyle \mathbf {a} (\mathbf {x} )}
is independent of
x
{\displaystyle \mathbf {x} }
, or in Fourier space,
a
^
(
k
)
=
0
{\displaystyle {\widehat {\mathbf {a} }}(\mathbf {k} )={\boldsymbol {0}}}
unless
k
=
0
{\displaystyle \mathbf {k} ={\boldsymbol {0}}}
.
The subspace
J
{\displaystyle {\mathcal {J}}}
of zero divergence, zero average value fields, i.e.,
∇
⋅
a
(
x
)
=
0
{\displaystyle {\boldsymbol {\nabla }}\cdot \mathbf {a} (\mathbf {x} )=0}
and
⟨
a
(
x
)
⟩
=
0
{\displaystyle \langle \mathbf {a} (\mathbf {x} )\rangle ={\boldsymbol {0}}}
, or in Fourier space,
a
^
(
0
)
=
0
{\displaystyle {\widehat {\mathbf {a} }}({\boldsymbol {0}})={\boldsymbol {0}}}
and
k
⋅
a
^
(
k
)
=
0
{\displaystyle \mathbf {k} \cdot {\widehat {\mathbf {a} }}(\mathbf {k} )=0}
.
The subspace
E
{\displaystyle {\mathcal {E}}}
of zero curl, zero average value fields, i.e.,
∇
×
a
(
x
)
=
0
{\displaystyle {\boldsymbol {\nabla }}\times \mathbf {a} (\mathbf {x} )={\boldsymbol {0}}}
and
⟨
a
(
x
)
⟩
=
0
{\displaystyle \langle \mathbf {a} (\mathbf {x} )\rangle ={\boldsymbol {0}}}
, or in Fourier space,
a
^
(
0
)
=
0
{\displaystyle {\widehat {\mathbf {a} }}({\boldsymbol {0}})={\boldsymbol {0}}}
and
a
^
(
k
)
=
α
(
k
)
k
{\displaystyle {\widehat {\mathbf {a} }}(\mathbf {k} )=\alpha (\mathbf {k} )~\mathbf {k} }
.
To determine the effective permittivity, we introduced the operator
Γ
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}}
with the properties
(1)
b
(
x
)
=
Γ
⋅
a
(
x
)
{\displaystyle {\text{(1)}}\qquad \mathbf {b} (\mathbf {x} )={\boldsymbol {\mathit {\Gamma }}}\cdot \mathbf {a} (\mathbf {x} )}
if and only if
b
∈
E
and
a
−
ϵ
0
⋅
b
∈
U
⊕
J
⟹
∇
⋅
(
a
−
ϵ
0
⋅
b
)
=
0
.
{\displaystyle \mathbf {b} \in {\mathcal {E}}\qquad {\text{and}}\qquad \mathbf {a} -{\boldsymbol {\epsilon }}_{0}\cdot \mathbf {b} \in {\mathcal {U}}\oplus {\mathcal {J}}\implies {\boldsymbol {\nabla }}\cdot (\mathbf {a} -{\boldsymbol {\epsilon }}_{0}\cdot \mathbf {b} )=0~.}
We also found that in Fourier space,
b
^
(
k
)
=
Γ
^
(
k
)
⋅
a
^
(
k
)
{\displaystyle {\widehat {\mathbf {b} }}(\mathbf {k} )={\boldsymbol {\mathit {\widehat {\Gamma }}}}(\mathbf {k} )\cdot {\widehat {\mathbf {a} }}(\mathbf {k} )}
where
(2)
Γ
^
(
k
)
=
{
k
⊗
k
k
⋅
ϵ
0
⋅
k
=
Γ
(
n
)
with
n
=
k
|
k
|
if
k
≠
0
0
if
k
=
0
.
{\displaystyle {\text{(2)}}\qquad {\boldsymbol {\mathit {\widehat {\Gamma }}}}(\mathbf {k} )={\begin{cases}{\cfrac {\mathbf {k} \otimes \mathbf {k} }{\mathbf {k} \cdot {\boldsymbol {\epsilon }}_{0}\cdot \mathbf {k} }}={\boldsymbol {\mathit {\Gamma }}}(\mathbf {n} )&\qquad {\text{with}}~~\mathbf {n} ={\cfrac {\mathbf {k} }{|\mathbf {k} |}}~~{\text{if}}~~\mathbf {k} \neq 0\\{\boldsymbol {\mathit {0}}}&\qquad {\text{if}}~~\mathbf {k} =0~.\end{cases}}}
Let us now derive a formula for the effective tensor. Recall that
the polarization is defined as
(3)
P
=
(
ϵ
−
ϵ
0
)
⋅
E
=
D
−
ϵ
0
⋅
E
{\displaystyle {\text{(3)}}\qquad \mathbf {P} =({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot \mathbf {E} =\mathbf {D} -{\boldsymbol {\epsilon }}_{0}\cdot \mathbf {E} }
where the permittivity tensor is being thought of as an operator that
operates on
E
{\displaystyle \mathbf {E} }
.
Also notice that
(4)
⟨
E
⟩
−
E
∈
E
(curl free)
.
{\displaystyle {\text{(4)}}\qquad \langle \mathbf {E} \rangle -\mathbf {E} \in {\mathcal {E}}\qquad {\text{(curl free)}}.}
From (3) we have
(5)
P
−
ϵ
0
⋅
[
⟨
E
⟩
−
E
]
=
D
−
ϵ
0
⋅
⟨
E
⟩
∈
U
⊕
J
(divergence free)
.
{\displaystyle {\text{(5)}}\qquad \mathbf {P} -{\boldsymbol {\epsilon }}_{0}\cdot [\langle \mathbf {E} \rangle -\mathbf {E} ]=\mathbf {D} -{\boldsymbol {\epsilon }}_{0}\cdot \langle \mathbf {E} \rangle \in {\mathcal {U}}\oplus {\mathcal {J}}\qquad {\text{(divergence free)}}.}
From the definition of
Γ
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}}
(equations (1) and
(2)) and using (4) and (5), we
can show that
(6)
Γ
⋅
P
=
⟨
E
⟩
−
E
.
{\displaystyle {\text{(6)}}\qquad {\boldsymbol {\mathit {\Gamma }}}\cdot \mathbf {P} =\langle \mathbf {E} \rangle -\mathbf {E} ~.}
Let
ϵ
−
ϵ
0
{\displaystyle {\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0}}
act on both sides of (6). Then
we get
(
ϵ
−
ϵ
0
)
⋅
Γ
⋅
P
=
(
ϵ
−
ϵ
0
)
⋅
⟨
E
⟩
−
(
ϵ
−
ϵ
0
)
⋅
E
=
(
ϵ
−
ϵ
0
)
⋅
⟨
E
⟩
−
P
.
{\displaystyle ({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot {\boldsymbol {\mathit {\Gamma }}}\cdot \mathbf {P} =({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot \langle \mathbf {E} \rangle -({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot \mathbf {E} =({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot \langle \mathbf {E} \rangle -\mathbf {P} ~.}
Therefore,
(7)
[
1
+
(
ϵ
−
ϵ
0
)
⋅
Γ
]
⋅
P
=
(
ϵ
−
ϵ
0
)
⋅
⟨
E
⟩
.
{\displaystyle {\text{(7)}}\qquad [{\boldsymbol {\mathit {1}}}+({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot {\boldsymbol {\mathit {\Gamma }}}]\cdot \mathbf {P} =({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot \langle \mathbf {E} \rangle ~.}
Inverting (7) gives
(8)
P
=
[
1
+
(
ϵ
−
ϵ
0
)
⋅
Γ
]
−
1
⋅
(
ϵ
−
ϵ
0
)
⋅
⟨
E
⟩
.
{\displaystyle {\text{(8)}}\qquad \mathbf {P} =[{\boldsymbol {\mathit {1}}}+({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot {\boldsymbol {\mathit {\Gamma }}}]^{-1}\cdot ({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot \langle \mathbf {E} \rangle ~.}
Averaging (8) gives
(9)
⟨
P
⟩
=
⟨
[
1
+
(
ϵ
−
ϵ
0
)
⋅
Γ
]
−
1
⋅
(
ϵ
−
ϵ
0
)
⟩
⋅
⟨
E
⟩
.
{\displaystyle {\text{(9)}}\qquad \langle \mathbf {P} \rangle =\langle [{\boldsymbol {\mathit {1}}}+({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot {\boldsymbol {\mathit {\Gamma }}}]^{-1}\cdot ({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\rangle \cdot \langle \mathbf {E} \rangle ~.}
Averaging (3) leads to the relation
(10)
⟨
P
⟩
=
(
ϵ
eff
−
ϵ
0
)
⋅
⟨
E
⟩
.
{\displaystyle {\text{(10)}}\qquad \langle \mathbf {P} \rangle =({\boldsymbol {\epsilon }}_{\text{eff}}-{\boldsymbol {\epsilon }}_{0})\cdot \langle \mathbf {E} \rangle ~.}
Comparing (9) and (10) shows us that
(11)
ϵ
eff
=
ϵ
0
+
⟨
[
1
+
(
ϵ
−
ϵ
0
)
⋅
Γ
]
−
1
⋅
(
ϵ
−
ϵ
0
)
⟩
.
{\displaystyle {\text{(11)}}\qquad {{\boldsymbol {\epsilon }}_{\text{eff}}={\boldsymbol {\epsilon }}_{0}+\langle [{\boldsymbol {\mathit {1}}}+({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot {\boldsymbol {\mathit {\Gamma }}}]^{-1}\cdot ({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\rangle ~.}}
Recall from the previous lecture that, in equation (11), the operator
ϵ
{\displaystyle {\boldsymbol {\epsilon }}}
is local in real space while the operator
Γ
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}}
is local in Fourier space. It is therefore not obvious how one can invert
[
1
+
(
ϵ
−
ϵ
0
)
⋅
Γ
]
{\displaystyle [{\boldsymbol {\mathit {1}}}+({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot {\boldsymbol {\mathit {\Gamma }}}]}
.
Let us define
δ
:=
ϵ
−
ϵ
0
.
{\displaystyle {\boldsymbol {\delta }}:={\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0}~.}
Then (8) can be written as
P
=
[
1
+
δ
⋅
Γ
]
−
1
⋅
δ
⋅
⟨
E
⟩
.
{\displaystyle \mathbf {P} =[{\boldsymbol {\mathit {1}}}+{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}}]^{-1}\cdot {\boldsymbol {\delta }}\cdot \langle \mathbf {E} \rangle ~.}
Assuming that
−
δ
⋅
Γ
<
1
{\displaystyle -{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}}<{\boldsymbol {\mathit {1}}}}
, we can expand the first operator in terms of an infinite series, i.e.,
[
1
+
δ
⋅
Γ
]
−
1
=
∑
n
=
0
∞
(
−
δ
⋅
Γ
)
n
.
{\displaystyle [{\boldsymbol {\mathit {1}}}+{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}}]^{-1}=\sum _{n=0}^{\infty }(-{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}})^{n}~.}
Then we have
P
=
[
∑
n
=
0
∞
(
−
δ
⋅
Γ
)
n
]
⋅
δ
⋅
⟨
E
⟩
.
{\displaystyle \mathbf {P} =\left[\sum _{n=0}^{\infty }(-{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}})^{n}\right]\cdot {\boldsymbol {\delta }}\cdot \langle \mathbf {E} \rangle ~.}
Also, from the definition of
P
{\displaystyle \mathbf {P} }
, we have
P
=
δ
⋅
E
.
{\displaystyle \mathbf {P} ={\boldsymbol {\delta }}\cdot \mathbf {E} ~.}
Hence,
δ
⋅
E
=
[
∑
n
=
0
∞
(
−
δ
⋅
Γ
)
n
]
⋅
δ
⋅
⟨
E
⟩
.
{\displaystyle {\boldsymbol {\delta }}\cdot \mathbf {E} =\left[\sum _{n=0}^{\infty }(-{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}})^{n}\right]\cdot {\boldsymbol {\delta }}\cdot \langle \mathbf {E} \rangle ~.}
Now,
(
−
δ
⋅
Γ
)
n
⋅
δ
=
(
−
1
)
n
δ
⋅
Γ
⋅
δ
⋅
Γ
…
δ
⋅
Γ
⋅
δ
=
δ
⋅
(
−
Γ
⋅
δ
)
n
.
{\displaystyle (-{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}})^{n}\cdot {\boldsymbol {\delta }}=(-1)^{n}~{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}}\cdot {\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}}\dots {\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}}\cdot {\boldsymbol {\delta }}={\boldsymbol {\delta }}\cdot (-{\boldsymbol {\mathit {\Gamma }}}\cdot {\boldsymbol {\delta }})^{n}~.}
Therefore,
δ
⋅
E
=
δ
⋅
[
∑
n
=
0
∞
(
−
Γ
⋅
δ
)
n
]
⋅
⟨
E
⟩
{\displaystyle {\boldsymbol {\delta }}\cdot \mathbf {E} ={\boldsymbol {\delta }}\cdot \left[\sum _{n=0}^{\infty }(-{\boldsymbol {\mathit {\Gamma }}}\cdot {\boldsymbol {\delta }})^{n}\right]\cdot \langle \mathbf {E} \rangle }
or,
(12)
E
=
[
∑
n
=
0
∞
(
−
Γ
⋅
δ
)
n
]
⋅
⟨
E
⟩
.
{\displaystyle {\text{(12)}}\qquad \mathbf {E} =\left[\sum _{n=0}^{\infty }(-{\boldsymbol {\mathit {\Gamma }}}\cdot {\boldsymbol {\delta }})^{n}\right]\cdot \langle \mathbf {E} \rangle ~.}
Let us now define
(13)
P
(
m
)
:=
[
∑
n
=
0
m
(
−
δ
⋅
Γ
)
n
]
⋅
δ
⋅
⟨
E
⟩
E
(
m
)
:=
[
∑
n
=
0
m
(
−
Γ
⋅
δ
)
n
]
⋅
⟨
E
⟩
.
{\displaystyle {\text{(13)}}\qquad {\begin{aligned}\mathbf {P} ^{(m)}&:=\left[\sum _{n=0}^{m}(-{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}})^{n}\right]\cdot {\boldsymbol {\delta }}\cdot \langle \mathbf {E} \rangle \\\mathbf {E} ^{(m)}&:=\left[\sum _{n=0}^{m}(-{\boldsymbol {\mathit {\Gamma }}}\cdot {\boldsymbol {\delta }})^{n}\right]\cdot \langle \mathbf {E} \rangle ~.\end{aligned}}}
Then we can write
(14)
Γ
⋅
P
(
m
)
=
⟨
E
⟩
−
E
(
m
+
1
)
and
δ
⋅
E
(
m
)
=
P
(
m
)
.
{\displaystyle {\text{(14)}}\qquad {\boldsymbol {\mathit {\Gamma }}}\cdot \mathbf {P} ^{(m)}=\langle \mathbf {E} \rangle -\mathbf {E} ^{(m+1)}\qquad {\text{and}}\qquad {\boldsymbol {\delta }}\cdot \mathbf {E} ^{(m)}=\mathbf {P} ^{(m)}~.}
These recurrence relations may be used to compute these fields inductively.
An algorithm that may be used is outlined below:
Set
m
=
0
{\displaystyle m=0}
. Then
E
(
0
)
=
⟨
E
⟩
.
{\displaystyle \mathbf {E} ^{(0)}=\langle \mathbf {E} \rangle ~.}
Compute
P
(
m
)
(
x
)
{\displaystyle \mathbf {P} ^{(m)}(\mathbf {x} )}
in real space using the relation
P
(
m
)
=
δ
⋅
E
(
m
)
.
{\displaystyle \mathbf {P} ^{(m)}={\boldsymbol {\delta }}\cdot \mathbf {E} ^{(m)}~.}
Take a fast Fourier transform to find
P
^
(
m
)
(
k
)
{\displaystyle {\widehat {\mathbf {P} }}^{(m)}(\mathbf {k} )}
.
From (14) we get
E
^
(
m
+
1
)
=
{
−
Γ
^
(
k
)
⋅
P
^
(
m
)
(
k
)
for
k
≠
0
⟨
E
⟩
for
k
=
0
.
{\displaystyle {\widehat {\mathbf {E} }}^{(m+1)}={\begin{cases}-{\boldsymbol {\mathit {\widehat {\Gamma }}}}(\mathbf {k} )\cdot {\widehat {\mathbf {P} }}^{(m)}(\mathbf {k} )&{\text{for}}~~\mathbf {k} \neq 0\\\langle \mathbf {E} \rangle &{\text{for}}~~\mathbf {k} =0~.\end{cases}}}
Compute
E
^
(
m
+
1
)
(
k
)
{\displaystyle {\widehat {\mathbf {E} }}^{(m+1)}(\mathbf {k} )}
in Fourier space.
Take an inverse fast Fourier transform to find
E
(
m
+
1
)
(
x
)
{\displaystyle \mathbf {E} ^{(m+1)}(\mathbf {x} )}
in real space.
Increment
m
{\displaystyle 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.
For simplicity, let us assume that
ϵ
0
{\displaystyle {\boldsymbol {\epsilon }}_{0}}
is isotropic, i.e.,
ϵ
0
=
ϵ
0
1
{\displaystyle {\boldsymbol {\epsilon }}_{0}=\epsilon _{0}~{\boldsymbol {\mathit {1}}}}
. Then,
Γ
=
Γ
1
ϵ
0
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}={\cfrac {{\boldsymbol {\mathit {\Gamma }}}_{1}}{\epsilon _{0}}}}
where
Γ
1
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{1}}
is the projection from
H
{\displaystyle {\mathcal {H}}}
onto
E
{\displaystyle {\mathcal {E}}}
.
Define the norm of a field
P
{\displaystyle \mathbf {P} }
as
|
P
|
:=
(
P
,
P
)
1
/
2
=
⟨
P
⟩
⋅
P
¯
1
/
2
.
{\displaystyle |\mathbf {P} |:=(\mathbf {P} ,\mathbf {P} )^{1/2}=\langle {\overline {\mathbf {P} \rangle \cdot \mathbf {P} }}^{1/2}~.}
Also, define the norm of a linear operator
A
{\displaystyle {\boldsymbol {A}}}
acting on
P
{\displaystyle \mathbf {P} }
as
‖
A
‖
:=
sup
P
∈
H
,
|
P
|
=
1
|
A
⋅
P
|
.
{\displaystyle \lVert {\boldsymbol {A}}\rVert _{}:=\sup _{\mathbf {P} \in {\mathcal {H}}~,~|\mathbf {P} |=1}|{\boldsymbol {A}}\cdot \mathbf {P} |~.}
Therefore,
|
A
⋅
P
|
≤
‖
A
‖
⋅
|
P
|
.
{\displaystyle |{\boldsymbol {A}}\cdot \mathbf {P} |\leq \lVert {\boldsymbol {A}}\rVert _{}\cdot |\mathbf {P} |~.}
Hence,
|
A
⋅
B
⋅
P
|
≤
‖
A
‖
⋅
|
B
⋅
P
|
≤
‖
A
‖
⋅
‖
B
‖
⋅
|
P
|
.
{\displaystyle |{\boldsymbol {A}}\cdot {\boldsymbol {B}}\cdot \mathbf {P} |\leq \lVert {\boldsymbol {A}}\rVert _{}\cdot |{\boldsymbol {B}}\cdot \mathbf {P} |\leq \lVert {\boldsymbol {A}}\rVert _{}\cdot \lVert {\boldsymbol {B}}\rVert \cdot |\mathbf {P} |~.}
So
‖
A
⋅
B
‖
≤
‖
A
‖
⋅
‖
B
‖
.
{\displaystyle \lVert {\boldsymbol {A}}\cdot {\boldsymbol {B}}\rVert _{}\leq \lVert {\boldsymbol {A}}\rVert _{}\cdot \lVert {\boldsymbol {B}}\rVert _{}~.}
In addition, we have the triangle inequality
‖
A
+
B
‖
≤
‖
A
‖
+
‖
B
‖
.
{\displaystyle \lVert {\boldsymbol {A}}+{\boldsymbol {B}}\rVert _{}\leq \lVert {\boldsymbol {A}}\rVert _{}+\lVert {\boldsymbol {B}}\rVert _{}~.}
So from (12) and (13), we have
|
E
−
E
(
m
)
|
=
|
∑
n
=
m
+
1
∞
(
−
δ
⋅
Γ
)
n
⋅
⟨
E
⟩
|
≤
‖
∑
n
=
m
+
1
∞
(
−
δ
⋅
Γ
)
n
‖
⋅
|
⟨
E
⟩
|
≤
(
∑
n
=
m
+
1
∞
‖
δ
ϵ
0
‖
n
⋅
‖
Γ
1
‖
n
)
⋅
|
⟨
E
⟩
|
.
{\displaystyle {\begin{aligned}|\mathbf {E} -\mathbf {E} ^{(m)}|&=\left|\sum _{n=m+1}^{\infty }(-{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}})^{n}\cdot \langle \mathbf {E} \rangle \right|\\&\leq \lVert \sum _{n=m+1}^{\infty }(-{\boldsymbol {\delta }}\cdot {\boldsymbol {\mathit {\Gamma }}})^{n}\rVert _{}\cdot |\langle \mathbf {E} \rangle |\\&\leq \left(\sum _{n=m+1}^{\infty }\lVert {\cfrac {\boldsymbol {\delta }}{\epsilon _{0}}}\rVert ^{n}\cdot \lVert {\boldsymbol {\mathit {\Gamma }}}_{1}\rVert _{}^{n}\right)\cdot |\langle \mathbf {E} \rangle |~.\end{aligned}}}
But
‖
Γ
1
‖
=
1
{\displaystyle \lVert {\boldsymbol {\mathit {\Gamma }}}_{1}\rVert _{}=1}
since it is a projection onto
E
{\displaystyle {\mathcal {E}}}
. Hence,
|
E
−
E
(
m
)
|
≤
(
∑
n
=
m
+
1
∞
‖
δ
ϵ
0
‖
n
)
⋅
|
⟨
E
⟩
|
.
{\displaystyle |\mathbf {E} -\mathbf {E} ^{(m)}|\leq \left(\sum _{n=m+1}^{\infty }\lVert {\cfrac {\boldsymbol {\delta }}{\epsilon _{0}}}\rVert ^{n}\right)\cdot |\langle \mathbf {E} \rangle |~.}
Therefore the series converges provided
‖
δ
ϵ
0
‖
<
1
.
{\displaystyle \lVert {\cfrac {\boldsymbol {\delta }}{\epsilon _{0}}}\rVert <1~.}
In that case
|
E
−
E
(
m
)
|
≤
γ
m
+
1
1
−
γ
|
⟨
E
⟩
|
{\displaystyle |\mathbf {E} -\mathbf {E} ^{(m)}|\leq {\cfrac {\gamma ^{m+1}}{1-\gamma }}~|\langle \mathbf {E} \rangle |}
where
γ
:=
‖
δ
ϵ
0
‖
.
{\displaystyle \gamma :=\lVert {\cfrac {\boldsymbol {\delta }}{\epsilon _{0}}}\rVert _{}~.}
To get a better understanding of the norm
γ
{\displaystyle \gamma }
, let us consider a
n
{\displaystyle n}
-phase composite with isotropic phases, i.e.,
ϵ
=
∑
i
=
1
n
χ
i
ϵ
i
1
{\displaystyle {\boldsymbol {\epsilon }}=\sum _{i=1}^{n}\chi _{i}~\epsilon _{i}~{\boldsymbol {\mathit {1}}}}
where
χ
i
=
{
1
in phase
i
0
otherwise
.
{\displaystyle \chi _{i}={\begin{cases}1&{\text{in phase}}~~i\\0&{\text{otherwise}}.\end{cases}}}
In this case,
γ
=
‖
δ
ϵ
0
‖
=
sup
P
∈
H
,
|
P
|
=
1
∑
i
=
1
n
⟨
(
ϵ
i
−
ϵ
0
ϵ
0
P
i
)
¯
⋅
(
ϵ
i
−
ϵ
0
ϵ
0
P
i
)
⟩
1
/
2
{\displaystyle \gamma =\lVert {\cfrac {\boldsymbol {\delta }}{\epsilon _{0}}}\rVert =\sup _{\mathbf {P} \in {\mathcal {H}}~,~|\mathbf {P} |=1}\sum _{i=1}^{n}\left\langle {\overline {\left({\cfrac {\epsilon _{i}-\epsilon _{0}}{\epsilon _{0}}}~\mathbf {P} _{i}\right)}}\cdot \left({\cfrac {\epsilon _{i}-\epsilon _{0}}{\epsilon _{0}}}~\mathbf {P} _{i}\right)\right\rangle ^{1/2}}
where
P
i
=
χ
i
P
and
|
P
|
=
(
∑
i
P
i
¯
⋅
P
)
1
/
2
=
1
.
{\displaystyle \mathbf {P} _{i}=\chi _{i}~\mathbf {P} \quad {\text{and}}\quad |\mathbf {P} |=\left(\sum _{i}{\overline {\mathbf {P} _{i}}}\cdot \mathbf {P} \right)^{1/2}=1~.}
Hence,
γ
=
sup
P
∈
H
,
|
P
|
=
1
∑
i
=
1
n
|
ϵ
i
−
ϵ
0
ϵ
0
|
⟨
P
¯
i
⋅
P
i
⟩
1
/
2
.
{\displaystyle \gamma =\sup _{\mathbf {P} \in {\mathcal {H}}~,~|\mathbf {P} |=1}\sum _{i=1}^{n}\left|{\cfrac {\epsilon _{i}-\epsilon _{0}}{\epsilon _{0}}}\right|\langle {\overline {\mathbf {P} }}_{i}\cdot \mathbf {P} _{i}\rangle ^{1/2}~.}
Since,
⟨
P
¯
i
⋅
P
i
⟩
1
/
2
{\displaystyle \langle {\overline {\mathbf {P} }}_{i}\cdot \mathbf {P} _{i}\rangle ^{1/2}}
are weights, it makes sense to
put the weights where
ϵ
i
{\displaystyle \epsilon _{i}}
is maximum. Hence, we can write
γ
=
max
i
|
ϵ
i
−
ϵ
0
ϵ
0
|
=
max
i
|
ϵ
i
−
ϵ
0
|
|
ϵ
0
|
.
{\displaystyle \gamma =\max _{i}\left|{\cfrac {\epsilon _{i}-\epsilon _{0}}{\epsilon _{0}}}\right|=\max _{i}{\cfrac {|\epsilon _{i}-\epsilon _{0}|}{|\epsilon _{0}|}}~.}
For
γ
{\displaystyle \gamma }
to be less than 1, we therefore require that, for all
i
{\displaystyle i}
,
|
ϵ
i
−
ϵ
0
|
<
|
ϵ
0
|
.
{\displaystyle |\epsilon _{i}-\epsilon _{0}|<|\epsilon _{0}|~.}
A geometrical representation of this situation is shown in Figure 1.
Figure 1. Allowable values of
ϵ
i
{\displaystyle \epsilon _{i}}
for convergence of series expansion.
If the value of
ϵ
0
{\displaystyle \epsilon _{0}}
is sufficiently large, then we get convergence
if all the
ϵ
i
{\displaystyle \epsilon _{i}}
s lie in one half of the complex plane (shown by
the green line in the figure).
Similarly, we can expand
ϵ
eff
=
ϵ
0
+
⟨
[
1
+
(
ϵ
−
ϵ
0
)
⋅
Γ
]
−
1
⋅
(
ϵ
−
ϵ
0
)
⟩
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}={\boldsymbol {\epsilon }}_{0}+\langle [{\boldsymbol {\mathit {1}}}+({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\cdot {\boldsymbol {\mathit {\Gamma }}}]^{-1}\cdot ({\boldsymbol {\epsilon }}-{\boldsymbol {\epsilon }}_{0})\rangle }
in the form
(15)
ϵ
eff
=
ϵ
0
+
∑
j
=
0
∞
Γ
0
⋅
δ
⋅
(
−
Γ
⋅
δ
)
j
⋅
Γ
0
{\displaystyle {\text{(15)}}\qquad {\boldsymbol {\epsilon }}_{\text{eff}}={\boldsymbol {\epsilon }}_{0}+\sum _{j=0}^{\infty }{\boldsymbol {\mathit {\Gamma }}}_{0}\cdot {\boldsymbol {\delta }}\cdot (-{\boldsymbol {\mathit {\Gamma }}}\cdot {\boldsymbol {\delta }})^{j}\cdot {\boldsymbol {\mathit {\Gamma }}}_{0}}
where
Γ
0
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{0}}
is a projection onto
U
{\displaystyle {\mathcal {U}}}
, i.e.,
Γ
0
⋅
P
=
⟨
P
⟩
.
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{0}\cdot \mathbf {P} =\langle \mathbf {P} \rangle ~.}
In this case, we find that the series converges provided
|
ϵ
i
−
ϵ
0
|
<
|
ϵ
0
|
for all
i
.
{\displaystyle |\epsilon _{i}-\epsilon _{0}|<|\epsilon _{0}|\qquad {\text{for all}}~~i~.}
Note that each term in (15) is an analytic function of
ϵ
1
{\displaystyle \epsilon _{1}}
(in fact a polynomial). So, if we truncate the series, we
have an analytic function of
ϵ
1
{\displaystyle \epsilon _{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
ϵ
eff
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}}
is an analytic
function of
ϵ
1
{\displaystyle \epsilon _{1}}
in the
ϵ
0
{\displaystyle \epsilon _{0}}
disk
(see Figure~1) with
r
<
|
ϵ
0
|
{\displaystyle r<|\epsilon _{0}|}
provided
|
ϵ
i
−
ϵ
0
|
<
ϵ
0
{\displaystyle |\epsilon _{i}-\epsilon _{0}|<\epsilon _{0}}
for
i
=
2
,
3
,
…
,
n
{\displaystyle i=2,3,\dots ,n}
.
Similarly, the effective tensor is an analytic function of
ϵ
2
{\displaystyle \epsilon _{2}}
,
ϵ
3
{\displaystyle \epsilon _{3}}
, etc.
Since
ϵ
eff
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}}
is independent of
ϵ
0
{\displaystyle \epsilon _{0}}
, by taking the union of
all such regions of analyticity, we conclude that
ϵ
eff
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}}
is an
analytic function of
ϵ
1
,
…
,
ϵ
n
{\displaystyle \epsilon _{1},\dots ,\epsilon _{n}}
provided all these
ϵ
i
{\displaystyle \epsilon _{i}}
s lie inside a half-plane (see Figure~1).
This means that there exists a
θ
{\displaystyle \theta }
such that
Re
(
e
i
θ
ϵ
j
)
>
0
for all
j
.
{\displaystyle {\text{Re}}(e^{i\theta }~\epsilon _{j})>0\quad {\text{for all}}~j~.}
A corollary of the above observations is the following. If each
ϵ
i
(
ω
)
{\displaystyle \epsilon _{i}(\omega )}
is an analytic function of
ω
{\displaystyle \omega }
for
Im
(
ω
)
>
0
{\displaystyle {\text{Im}}(\omega )>0}
(which is what one expects with
ω
{\displaystyle \omega }
as the frequency) and
Im
[
ϵ
i
(
ω
)
]
>
0
{\displaystyle {\text{Im}}[\epsilon _{i}(\omega )]>0}
for all
ω
{\displaystyle \omega }
with
Re
ω
>
0
{\displaystyle {\text{Re}}{\omega }>0}
,
then
ϵ
eff
(
ω
)
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}(\omega )}
will be analytic in
ω
{\displaystyle \omega }
.
Now, if
D
=
ϵ
⋅
E
{\displaystyle \mathbf {D} ={\boldsymbol {\epsilon }}\cdot \mathbf {E} }
, we have
λ
D
=
λ
ϵ
⋅
E
for all constants
λ
.
{\displaystyle \lambda ~\mathbf {D} =\lambda ~{\boldsymbol {\epsilon }}\cdot \mathbf {E} \qquad {\text{for all constants}}~\lambda ~.}
Therefore,
⟨
D
⟩
=
ϵ
eff
⋅
⟨
E
⟩
⟹
⟨
λ
D
⟩
=
λ
ϵ
eff
⋅
⟨
E
⟩
.
{\displaystyle \langle \mathbf {D} \rangle ={\boldsymbol {\epsilon }}_{\text{eff}}\cdot \langle \mathbf {E} \rangle \qquad \implies \qquad \langle \lambda ~\mathbf {D} \rangle =\lambda ~{\boldsymbol {\epsilon }}_{\text{eff}}\cdot \langle \mathbf {E} \rangle ~.}
This means that
(
λ
ϵ
)
eff
=
λ
ϵ
eff
.
{\displaystyle (\lambda ~{\boldsymbol {\epsilon }})_{\text{eff}}=\lambda ~{\boldsymbol {\epsilon }}_{\text{eff}}~.}
Therefore,
ϵ
eff
(
ϵ
1
,
ϵ
2
,
…
,
ϵ
n
)
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}(\epsilon _{1},\epsilon _{2},\dots ,\epsilon _{n})}
is homogeneous of degree one and
ϵ
eff
(
λ
ϵ
1
,
λ
ϵ
2
,
…
,
λ
ϵ
n
)
=
λ
ϵ
eff
(
ϵ
1
,
ϵ
2
,
…
,
ϵ
n
)
.
{\displaystyle {{\boldsymbol {\epsilon }}_{\text{eff}}(\lambda ~\epsilon _{1},\lambda ~\epsilon _{2},\dots ,\lambda ~\epsilon _{n})=\lambda ~{\boldsymbol {\epsilon }}_{\text{eff}}(\epsilon _{1},\epsilon _{2},\dots ,\epsilon _{n})~.}}
For a two-phase composite, if we take
λ
=
1
ϵ
2
{\displaystyle \lambda ={\cfrac {1}{\epsilon _{2}}}}
we get
ϵ
eff
(
ϵ
1
,
ϵ
2
)
=
ϵ
2
ϵ
eff
(
ϵ
1
/
ϵ
2
,
1
)
.
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}(\epsilon _{1},\epsilon _{2})=\epsilon _{2}~{\boldsymbol {\epsilon }}_{\text{eff}}(\epsilon _{1}/\epsilon _{2},1)~.}
Therefore, it suffices to study the analytic function
ϵ
eff
(
ϵ
1
,
1
)
=
ϵ
eff
(
ϵ
1
)
{\displaystyle {\boldsymbol {\epsilon }}_{\text{eff}}(\epsilon _{1},1)={\boldsymbol {\epsilon }}_{\text{eff}}(\epsilon _{1})}
. For further
details see Milton02 .
[Milton02] G. W. Milton. Theory of Composites . Cambridge University Press, New York, 2002.
[Mouli94] 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.
[Mouli98] 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.
[Rudin76] W. Rudin. Principles of Mathematical Analysis . McGraw-Hill, New York, 1976.