text:
Waves in composites and metamaterials/Bloch waves in elastodynamics and bubbly fluids

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 discussed Bloch wave solutions for Maxwell's equations in periodic media in the quasistatic limit. Recall the periodic medium shown in Figure 1. The lattice spacing is $$\eta$$.

The permittivity and permeability of the medium are periodic functions of the form

\boldsymbol{\epsilon}(\mathbf{y} + \mathbf{a}_i) = \boldsymbol{\epsilon}(\mathbf{y}) ~; \boldsymbol{\mu}(\mathbf{y} + \mathbf{a}_i) = \boldsymbol{\mu}(\mathbf{y}) $$ where $$\mathbf{a}_i$$ are the primitive lattice vectors. We defined periodic functions

\begin{align} \boldsymbol{\epsilon}_\eta(\mathbf{x}) &:= \boldsymbol{\epsilon}\left(\cfrac{\mathbf{x}}{\eta}\right) ~; \boldsymbol{\mu}_\eta(\mathbf{x}) := \boldsymbol{\mu}\left(\cfrac{\mathbf{x}}{\eta}\right) ~; \\ \mathbf{D}_\eta(\mathbf{x}) &:= \mathbf{D}\left(\cfrac{\mathbf{x}}{\eta}\right) ~; \mathbf{B}_\eta(\mathbf{x}) := \mathbf{B}\left(\cfrac{\mathbf{x}}{\eta}\right) ~; \mathbf{E}_\eta(\mathbf{x}) := \mathbf{E}\left(\cfrac{\mathbf{x}}{\eta}\right) ~; \mathbf{H}_\eta(\mathbf{x}) := \mathbf{H}\left(\cfrac{\mathbf{x}}{\eta}\right) ~. \end{align} $$ We then wrote Maxwell's equations (at fixed frequency) as
 * $$ \text{(1)} \qquad

{ \boldsymbol{\nabla} \cdot \mathbf{D}_\eta = 0 ~; \boldsymbol{\nabla} \cdot \mathbf{B}_\eta = 0 ~; \boldsymbol{\nabla} \times \mathbf{E}_\eta - i\omega~\mathbf{B}_\eta = \boldsymbol{0} ~; \boldsymbol{\nabla} \times \mathbf{H}_\eta + i\omega~\mathbf{D}_\eta = \boldsymbol{0} ~. } $$ We also found that the constitutive relations could be expressed as
 * $$ \text{(2)} \qquad

{ \mathbf{d}_\eta(\mathbf{x}) = \boldsymbol{\epsilon}_\eta(\mathbf{x})\cdot\mathbf{e}_\eta(\mathbf{x}) ~; \mathbf{b}_\eta(\mathbf{x}) = \boldsymbol{\mu}_\eta(\mathbf{x})\cdot\mathbf{h}_\eta(\mathbf{x}) ~. } $$ Next, we examined solutions to (1) that correspond to Bloch waves with wavevector $$\mathbf{k}$$ (possibly complex), i.e.,

\mathbf{E}_\eta(\mathbf{x}) = e^{i\mathbf{k}\cdot\mathbf{x}}~\mathbf{e}_\eta(\mathbf{x}) ~; \mathbf{D}_\eta(\mathbf{x}) = e^{i\mathbf{k}\cdot\mathbf{x}}~\mathbf{d}_\eta(\mathbf{x}) ~; \mathbf{H}_\eta(\mathbf{x}) = e^{i\mathbf{k}\cdot\mathbf{x}}~\mathbf{h}_\eta(\mathbf{x}) ~; \mathbf{B}_\eta(\mathbf{x}) = e^{i\mathbf{k}\cdot\mathbf{x}}~\mathbf{b}_\eta(\mathbf{x}) ~. $$ The wavelength of the Bloch waves is

\lambda = \cfrac{2\pi}{|\text{Re}(\mathbf{k})|} ~. $$ For the quasistatic limit, we look for solutions that satisfy

\cfrac{2\pi}{|\mathbf{k}|} \gg |\eta~\mathbf{a}_i| ~. $$ Such a condition is satisfied if $$\eta\rightarrow 0$$. Following standard multiple scale analysis, we assumed that the periodic complex fields have the expansions
 * $$ \text{(3)} \qquad

\begin{align} \mathbf{e}_\eta(\mathbf{x}) & = \mathbf{e}_0(\mathbf{y}) + \eta~\mathbf{e}_1(\mathbf{y}) + \eta^2~\mathbf{e}_2(\mathbf{y}) + \dots \\ \mathbf{d}_\eta(\mathbf{x}) & = \mathbf{d}_0(\mathbf{y}) + \eta~\mathbf{d}_1(\mathbf{y}) + \eta^2~\mathbf{d}_2(\mathbf{y}) + \dots \\ \mathbf{h}_\eta(\mathbf{x}) & = \mathbf{h}_0(\mathbf{y}) + \eta~\mathbf{h}_1(\mathbf{y}) + \eta^2~\mathbf{h}_2(\mathbf{y}) + \dots \\ \mathbf{b}_\eta(\mathbf{x}) & = \mathbf{b}_0(\mathbf{y}) + \eta~\mathbf{b}_1(\mathbf{y}) + \eta^2~\mathbf{b}_2(\mathbf{y}) + \dots \end{align} $$ and that dependence of $$\omega$$ on $$\eta$$ and $$\mathbf{k}$$ has an expansion of the form
 * $$ \text{(4)} \qquad

\omega = \omega_\eta^j(\mathbf{k}) = \omega_0^j(\mathbf{y}) + \eta~\omega_1^j(\mathbf{y}) + \eta^2~\omega_2^j(\mathbf{y}) + \dots $$ Defining

\boldsymbol{\nabla}_y := \left(\frac{\partial }{\partial y_1}, \frac{\partial }{\partial y_2}, \frac{\partial }{\partial y_3}\right) $$ we then showed that we got the equations in the quasistatic limit
 * $$ \text{(5)} \qquad

{ \boldsymbol{\nabla}_y \cdot\mathbf{d}_0(\mathbf{y}) = 0 ~; \boldsymbol{\nabla}_y \cdot\mathbf{b}_0(\mathbf{y}) = 0 ~; \boldsymbol{\nabla}_y\times\mathbf{e}_0(\mathbf{y}) = \boldsymbol{0} ~; \boldsymbol{\nabla}_y\times\mathbf{h}_0(\mathbf{y}) = \boldsymbol{0} } $$ and the constitutive equations
 * $$ \text{(6)} \qquad

{ \mathbf{d}_0(\mathbf{y}) = \boldsymbol{\epsilon}(\mathbf{y})\cdot\mathbf{e}_0(\mathbf{y}) ~; \mathbf{b}_0(\mathbf{y}) = \boldsymbol{\mu}(\mathbf{y})\cdot\mathbf{h}_0(\mathbf{y}) ~. } $$ Note that the quasistatic equations (5), (6) have the same form as the Maxwell equations (1) and (2).

We also found that the volume averaged fields ($$\langle \bullet \rangle$$) satisfy the equations
 * $$ \text{(7)} \qquad

{   \langle \boldsymbol{\nabla}_y \cdot\mathbf{d}_1(\mathbf{y})\rangle = 0 ~; \langle \boldsymbol{\nabla}_y \cdot\mathbf{b}_1(\mathbf{y})\rangle = 0 ~; \langle \boldsymbol{\nabla}_y\times\mathbf{e}_1(\mathbf{y})\rangle = \boldsymbol{0} ~; \langle \boldsymbol{\nabla}_y\times\mathbf{h}_1(\mathbf{y})\rangle = \boldsymbol{0} ~. } $$ along with the necessary conditions
 * $$ \text{(8)} \qquad

{   \mathbf{k}\times\langle \mathbf{e}_0(\mathbf{y}) \rangle - \omega_0^j~\langle \mathbf{b}_0(\mathbf{y}) \rangle = \boldsymbol{0} ~; \mathbf{k}\times\langle \mathbf{h}_0(\mathbf{y}) \rangle + \omega_0^j~\langle \mathbf{d}_0(\mathbf{y}) \rangle = \boldsymbol{0} ~. } $$

Effective Permittivity and Permeability
Let $$\boldsymbol{\epsilon}_{\text{eff}}$$ be an effective permittivity tensor associated with the field $$\boldsymbol{\epsilon}(\mathbf{y})$$ and let $$\boldsymbol{\mu}_{\text{eff}}$$ be an effective permeability tensor associated with the field $$\boldsymbol{\mu}(\mathbf{y})$$. These effective tensors are defined through
 * $$ \text{(9)} \qquad

\langle \mathbf{d}_0 \rangle = \boldsymbol{\epsilon}_{\text{eff}} \cdot \langle \mathbf{e}_0 \rangle ~; \langle \mathbf{b}_0 \rangle = \boldsymbol{\mu}_{\text{eff}} \cdot \langle \mathbf{h}_0 \rangle ~. $$ Plugging (9) into (8) gives
 * $$ \text{(10)} \qquad

\mathbf{k}\times\langle \mathbf{e}_0 \rangle - \omega_0^j~\boldsymbol{\mu}_{\text{eff}} \cdot \langle \mathbf{h}_0 \rangle = \boldsymbol{0} ~; \mathbf{k}\times\langle \mathbf{h}_0 \rangle + \omega_0^j~\boldsymbol{\epsilon}_{\text{eff}} \cdot \langle \mathbf{e}_0 \rangle = \boldsymbol{0} ~. $$ Eliminate $$\langle \mathbf{h}_0 \rangle$$ from (10) to get

\cfrac{1}{\omega_0^j}~ \mathbf{k}\times\left\{\boldsymbol{\mu}_{\text{eff}}^{-1}\cdot[\mathbf{k}\times\langle \mathbf{e}_0 \rangle]\right\} + \omega_0^j~\boldsymbol{\epsilon}_{\text{eff}} \cdot \langle \mathbf{e}_0 \rangle = \boldsymbol{0} $$ or,
 * $$ \text{(11)} \qquad

\boldsymbol{\epsilon}_{\text{eff}}^{-1}\cdot\left( \mathbf{k}\times\left\{\boldsymbol{\mu}_{\text{eff}}^{-1}\cdot[\mathbf{k}\times\langle \mathbf{e}_0 \rangle]\right\}\right) + (\omega_0^j)^2~\langle \mathbf{e}_0 \rangle = \boldsymbol{0} ~. $$ Define the matrix $$\boldsymbol{A}_{\text{eff}}(\mathbf{k})$$ such that for all vectors $$\mathbf{v}$$
 * $$ \text{(12)} \qquad

\boldsymbol{A}_{\text{eff}}(\mathbf{k})\cdot\mathbf{v} = \boldsymbol{\epsilon}_{\text{eff}}^{-1}\cdot\left( \mathbf{k}\times\left\{\boldsymbol{\mu}_{\text{eff}}^{-1}\cdot[\mathbf{k}\times\mathbf{v}]\right\}\right) ~. $$ Then (11) can be written as
 * $$ \text{(13)} \qquad

{ \boldsymbol{A}_{\text{eff}}(\mathbf{k})\cdot\langle \mathbf{e}_0 \rangle + (\omega_0^j)^2~\langle \mathbf{e}_0 \rangle = \boldsymbol{0} ~. } $$ The eigenvalue problem (13) may be used to determine the dispersion relations $$\omega_0^j$$ and also possible values of $$\langle \mathbf{e}_0 \rangle$$ associated with each Block wave mode. Once we know $$\langle \mathbf{e}_0 \rangle$$ we can then find $$\langle \mathbf{h}_0 \rangle$$.

Note that the eigenvalue of $$\boldsymbol{A}_{\text{eff}}$$ is zero when $$\mathbf{e}_0 = \mathbf{k}$$. So, in practice, it is necessary to examine only the other two eigenvalues.

 Summary:

If the average field $$\langle \mathbf{e}_0 \rangle$$ is to be a Bloch wave solution of wavevector $$\mathbf{k}$$ and frequency $$\omega_0^j$$, then $$\langle \mathbf{e}_0 \rangle$$ must be an eigenvector of the matrix $$\boldsymbol{A}_{\text{eff}}(\mathbf{k})$$ and $$(\omega_0^j)^2$$ must be the associated eigenvalue. More general solutions can be obtained by superposing different Block wave solutions.

Isotropic constituents
When the permittivity and permeability of the medium are isotropic, equation (12) becomes

\boldsymbol{A}_{\text{eff}}(\mathbf{k})\cdot\mathbf{v} = \cfrac{1}{\epsilon_{\text{eff}}~\mu_{\text{eff}}}\cdot\left(\mathbf{k}\times\mathbf{k}\times\mathbf{v}\right) $$ or,
 * $$ \text{(14)} \qquad

\boldsymbol{A}_{\text{eff}}(\mathbf{k})\cdot\mathbf{v} = \cfrac{1}{\epsilon_{\text{eff}}~\mu_{\text{eff}}}\cdot\left[(\mathbf{k}\cdot\mathbf{k})~\mathbf{v} - (\mathbf{k}\cdot\mathbf{v})~\mathbf{k} \right] ~. $$ Since $$\mathbf{k}\cdot\langle \mathbf{e}_0 \rangle = 0$$, equation (13) becomes
 * $$ \text{(15)} \qquad

\cfrac{1}{\epsilon_{\text{eff}}~\mu_{\text{eff}}}~(\mathbf{k}\cdot\mathbf{k}) + (\omega_0^j)^2 = 0 ~. $$ If we assume that

\mathbf{k} = \left(\cfrac{2\pi}{\lambda} + \cfrac{i}{\delta}\right)~\mathbf{n} $$ where $$\lambda$$ is the wavelength, $$\delta$$ is the attenuation length, and $$\mathbf{n}$$ is a real unit vector, equation (15) gives

\cfrac{1}{\epsilon_{\text{eff}}~\mu_{\text{eff}}}~ \left(\cfrac{2\pi}{\lambda} + \cfrac{i}{\delta}\right)^2 + (\omega_0^j)^2 = 0 ~. $$ Therefore for a given frequency, the wavelength and attenuation length are determined by the complex permittivity $$\epsilon_{\text{eff}}$$ and the complex permeability $$\mu_{\text{eff}}$$.

Elastic wave propagation in the quasistatic limit
Let us now perform the exercise for linear elastic materials. Recall that the momentum equation at fixed frequency is given by

\boldsymbol{\nabla} \cdot \boldsymbol{\sigma}(\mathbf{x}) = - \omega^2~\rho(\mathbf{x})~\mathbf{u}(\mathbf{x}) $$ where $$\boldsymbol{\sigma}$$ is the stress, $$\rho$$ is the mass density, and $$\mathbf{u}$$ is the displacement. For infinitesimal deformations, the strain is related to the displacement by

\boldsymbol{\epsilon}(\mathbf{x}) = \frac{1}{2}\left\{\boldsymbol{\nabla}\mathbf{u}(\mathbf{x}) + [\boldsymbol{\nabla}\mathbf{u}(\mathbf{x})]^T\right\} ~. $$ Also, the stress is related to the strain for linear elastic materials by

\boldsymbol{\sigma}(\mathbf{x}) = \boldsymbol{\mathsf{C}}(\mathbf{x}):\boldsymbol{\epsilon}(\mathbf{x}) ~. $$ For a periodic medium, define

\boldsymbol{\sigma}_\eta(\mathbf{x}) = \boldsymbol{\sigma}(\mathbf{x}/\eta) ~; \boldsymbol{\epsilon}_\eta(\mathbf{x}) = \boldsymbol{\epsilon}(\mathbf{x}/\eta) ~; \mathbf{u}_\eta(\mathbf{x}) = \mathbf{u}(\mathbf{x}/\eta) ~; \rho_\eta(\mathbf{x}) = \rho(\mathbf{x}/\eta) ~; \boldsymbol{\mathsf{C}}_\eta(\mathbf{x}) = \boldsymbol{\mathsf{C}}(\mathbf{x}/\eta) ~. $$ Then the governing equations can be written as
 * $$ \text{(16)} \qquad

\boldsymbol{\nabla} \cdot \boldsymbol{\sigma}_\eta(\mathbf{x}) = - \omega^2~\rho_\eta(\mathbf{x})~\mathbf{u}_\eta(\mathbf{x}) ~; \boldsymbol{\epsilon}_\eta(\mathbf{x}) = \frac{1}{2}\left\{\boldsymbol{\nabla}\mathbf{u}_\eta(\mathbf{x}) + [\boldsymbol{\nabla}\mathbf{u}_\eta(\mathbf{x})]^T\right\} ~; \boldsymbol{\sigma}_\eta(\mathbf{x}) = \boldsymbol{\mathsf{C}}_\eta(\mathbf{x}):\boldsymbol{\epsilon}_\eta(\mathbf{x}) ~. $$ Let us examine Bloch wave solutions to these equations of the form
 * $$ \text{(17)} \qquad

\boldsymbol{\sigma}_\eta(\mathbf{x}) = e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\sigma}}_\eta(\mathbf{x}) ~; \boldsymbol{\epsilon}_\eta(\mathbf{x}) = e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\epsilon}}_\eta(\mathbf{x}) ~; \mathbf{u}_\eta(\mathbf{x}) = e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\mathbf{u}}_\eta(\mathbf{x}) ~. $$ Plug (17) into (16) to get

\begin{align} \boldsymbol{\nabla} \cdot \left(e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\sigma}}_\eta(\mathbf{x})\right) & = - \omega^2~\rho_\eta(\mathbf{x})~ e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\mathbf{u}}_\eta(\mathbf{x}) \\ e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\epsilon}}_\eta(\mathbf{x}) & = \frac{1}{2}\left\{ \boldsymbol{\nabla} \left(e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\mathbf{u}}_\eta(\mathbf{x})\right) + \left[\boldsymbol{\nabla} \left(e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\mathbf{u}}_\eta(\mathbf{x})\right)\right] \right\} \\ e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\sigma}}_\eta(\mathbf{x}) & = \boldsymbol{\mathsf{C}}_\eta(\mathbf{x}): \left(e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\epsilon}}_\eta(\mathbf{x}) \right) \end{align} $$ or,

\begin{align} & e^{i~\mathbf{k}\cdot\mathbf{x}}~\left(\boldsymbol{\nabla} \cdot \widehat{\boldsymbol{\sigma}}_\eta +    i~\widehat{\boldsymbol{\sigma}}_\eta~\cdot\mathbf{k}\right) = - e^{i~\mathbf{k}\cdot\mathbf{x}}~\omega^2~\rho_\eta~\widehat{\mathbf{u}}_\eta \\ & e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\epsilon}}_\eta = e^{i~\mathbf{k}\cdot\mathbf{x}}~\left[ \cfrac{i}{2}\left(\widehat{\mathbf{u}}_\eta\otimes\mathbf{k} + \mathbf{k}\otimes\widehat{\mathbf{u}}_\eta      \right) + \frac{1}{2}~\left(\boldsymbol{\nabla} \widehat{\mathbf{u}}_\eta + [\boldsymbol{\nabla} \widehat{\mathbf{u}}_\eta]^T\right) \right] \\ & e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\sigma}}_\eta = e^{i~\mathbf{k}\cdot\mathbf{x}}~\boldsymbol{\mathsf{C}}_\eta: \widehat{\boldsymbol{\epsilon}}_\eta \end{align} $$ or,
 * $$ \text{(18)} \qquad

{ \begin{align} &\boldsymbol{\nabla} \cdot \widehat{\boldsymbol{\sigma}}_\eta + i~\widehat{\boldsymbol{\sigma}}_\eta~\cdot\mathbf{k} + \omega^2~\rho_\eta~\widehat{\mathbf{u}}_\eta = 0 \\ &\widehat{\boldsymbol{\epsilon}}_\eta = \cfrac{i}{2}\left(\widehat{\mathbf{u}}_\eta\otimes\mathbf{k} + \mathbf{k}\otimes\widehat{\mathbf{u}}_\eta      \right) + \frac{1}{2}~\left(\boldsymbol{\nabla} \widehat{\mathbf{u}}_\eta + [\boldsymbol{\nabla} \widehat{\mathbf{u}}_\eta]^T\right) \\ &\widehat{\boldsymbol{\sigma}}_\eta = \boldsymbol{\mathsf{C}}_\eta: \widehat{\boldsymbol{\epsilon}}_\eta ~. \end{align} } $$ Equations (18) have solutions if the frequency $$\omega$$ takes one of a discrete set of values, $$\omega_\eta^j(\mathbf{k})$$, $$j= 1, 2, \dots$$. Let us assume that the fields have perturbation expansions in powers of $$\eta$$, i.e,
 * $$ \text{(19)} \qquad

\begin{align} &\widehat{\boldsymbol{\sigma}}_\eta(\mathbf{x}) = \boldsymbol{\sigma}_0(\mathbf{y}) + \eta~\boldsymbol{\sigma}_1(\mathbf{y}) + \eta^2~\boldsymbol{\sigma}_2(\mathbf{y}) + \dots \\ &\widehat{\boldsymbol{\epsilon}}_\eta(\mathbf{x}) = \boldsymbol{\epsilon}_0(\mathbf{y}) + \eta~\boldsymbol{\epsilon}_1(\mathbf{y}) + \eta^2~\boldsymbol{\epsilon}_2(\mathbf{y}) + \dots \\ &\widehat{\mathbf{u}}_\eta(\mathbf{x}) = \mathbf{u}_0(\mathbf{y}) + \eta~\mathbf{u}_1(\mathbf{y}) + \eta^2~\mathbf{u}_2(\mathbf{y}) + \dots \end{align} $$ where the functions $$\boldsymbol{\sigma}_i$$, $$\boldsymbol{\epsilon}_i$$, $$\mathbf{u}_i$$ are periodic. Let us also assume that

\omega = \omega_\eta^j(\mathbf{k}) = \omega_0^j(\mathbf{y}) + \eta~\omega_1^j(\mathbf{y}) + \eta^2~\omega_2^j(\mathbf{y}) + \dots $$ Plugging these expansions into (18) and using the definition

\boldsymbol{\nabla}_y := \left(\frac{\partial }{\partial y_1}, \frac{\partial }{\partial y_2}, \frac{\partial }{\partial y_3}\right) \qquad \implies \boldsymbol{\nabla} = \cfrac{1}{\eta}~\boldsymbol{\nabla}_y \quad \text{since} \quad \mathbf{y} = \cfrac{\mathbf{x}}{\eta} $$ we get
 * $$ \text{(20)} \qquad

\begin{align} &\cfrac{1}{\eta}\boldsymbol{\nabla}_y \cdot(\boldsymbol{\sigma}_0 + \eta~\boldsymbol{\sigma}_1 + \dots) + i~(\boldsymbol{\sigma}_0 + \eta~\boldsymbol{\sigma}_1 + \dots)~\cdot\mathbf{k} + \rho_\eta~ (\omega_0^j + \eta~\omega_1^j + \dots)^2~ (\mathbf{u}_0 + \eta~\mathbf{u}_1 + \dots) = 0 \\ & \boldsymbol{\epsilon}_0 + \eta~\boldsymbol{\epsilon}_1 + \dots = \cfrac{i}{2}\left[(\mathbf{u}_0 + \eta~\mathbf{u}_1 + \dots)\otimes\mathbf{k} + \mathbf{k}\otimes(\mathbf{u}_0 + \eta~\mathbf{u}_1 + \dots) \right] +\\ & \qquad\qquad\qquad\qquad \cfrac{1}{2\eta}~\left(\boldsymbol{\nabla}_y (\mathbf{u}_0 + \eta~\mathbf{u}_1 + \dots)+                [\boldsymbol{\nabla}_y (\mathbf{u}_0 + \eta~\mathbf{u}_1 + \dots)]^T\right) \\ & \boldsymbol{\sigma}_0 + \eta~\boldsymbol{\sigma}_1 + \dots = \boldsymbol{\mathsf{C}}_\eta: (\boldsymbol{\epsilon}_0 + \eta~\boldsymbol{\epsilon}_1 + \dots) ~. \end{align} $$ Collecting terms containing $$1/\eta$$ from (20) we get
 * $$\begin{align}

& {\boldsymbol{\nabla}_y \cdot\boldsymbol{\sigma}_0 = 0 } \text{(21)} \qquad \\ & {\boldsymbol{\nabla}_y \mathbf{u}_0 + [\boldsymbol{\nabla}_y \mathbf{u}_0]^T = 0}  ~. \text{(22)} \qquad \end{align}$$ Collecting terms of order $$\eta^0$$ from (20) we get
 * $$\begin{align}

& { \boldsymbol{\nabla}_y \cdot\boldsymbol{\sigma}_1 + i~\boldsymbol{\sigma}_0\cdot\mathbf{k} + \rho_\eta~(\omega_0^j)^2~\mathbf{u}_0 = 0} \text{(23)} \qquad \\ & { \boldsymbol{\epsilon}_0 = \cfrac{i}{2}~(\mathbf{u}_0\otimes\mathbf{k} + \mathbf{k}\otimes\mathbf{u}_0) + \frac{1}{2}~(\boldsymbol{\nabla}\mathbf{u}_1 + [\boldsymbol{\nabla}\mathbf{u}_1]^T)} \text{(24)} \qquad \\ & { \boldsymbol{\sigma}_0 = \boldsymbol{\mathsf{C}}_\eta:\boldsymbol{\epsilon}_0} \text{(25)} \qquad ~. \end{align}$$ From (22) we see that $$\mathbf{u}_0(\mathbf{y})$$ must be linear in $$\mathbf{y}$$. We also know that (from our definition), $$\mathbf{u}_0(\mathbf{y})$$ is periodic in $$\mathbf{y}$$. A function that is both linear and periodic must be constant. Therefore,

\mathbf{u}_0(\mathbf{y}) = \widehat{\mathbf{u}}_0 = ~\text{constant} ~. $$ Then we can write (24) in the form
 * $$ \text{(26)} \qquad

\boldsymbol{\epsilon}_0 = \cfrac{i}{2}~(\widehat{\mathbf{u}}_0\otimes\mathbf{k} + \mathbf{k}\otimes\widehat{\mathbf{u}}_0) + \frac{1}{2}~(\boldsymbol{\nabla}\mathbf{u}_1 + [\boldsymbol{\nabla}\mathbf{u}_1]^T) \qquad \implies { \boldsymbol{\epsilon}_0 = \frac{1}{2}~\left(\boldsymbol{\nabla}_y \widehat{\mathbf{u}}(\mathbf{y}) + [\boldsymbol{\nabla}_y \widehat{\mathbf{u}}(\mathbf{y})]^T\right) } $$ where

\widehat{\mathbf{u}}(\mathbf{y}) = i~(\mathbf{k}\cdot\mathbf{y})~\widehat{\mathbf{u}}_0 + \mathbf{u}_1(\mathbf{y}) ~. $$ Note that equations (21), (26), and (25) have a form similar to that of the elasticity equations in the absence of inertial and body forces.

The complex effective elasticity tensor may be defined via the relation
 * $$ \text{(27)} \qquad

{ \langle \boldsymbol{\sigma}_0 \rangle = \boldsymbol{\mathsf{C}}_\text{eff} : \langle \boldsymbol{\epsilon}_0 \rangle ~. } $$ Taking the average of (26) over the periodic cell, we get

\langle \boldsymbol{\epsilon}_0 \rangle = \cfrac{i}{2}~(\widehat{\mathbf{u}}_0\otimes\mathbf{k} + \mathbf{k}\otimes\widehat{\mathbf{u}}_0) + \frac{1}{2}~\left(\langle \boldsymbol{\nabla} \mathbf{u}_1\rangle + \langle [\boldsymbol{\nabla} \mathbf{u}_1]^T\rangle\right) $$ or
 * $$ \text{(28)} \qquad

{ \langle \boldsymbol{\epsilon}_0 \rangle = \cfrac{i}{2}~(\widehat{\mathbf{u}}_0\otimes\mathbf{k} + \mathbf{k}\otimes\widehat{\mathbf{u}}_0)~. } $$ Also, taking the average of (23) we get

\langle \boldsymbol{\nabla}_y \cdot\boldsymbol{\sigma}_1 \rangle + i~\langle \boldsymbol{\sigma}_0 \rangle\cdot\mathbf{k} + \langle \rho_\eta \rangle~(\omega_0^j)^2~\widehat{\mathbf{u}}_0 = 0 $$ or,
 * $$ \text{(29)} \qquad

{ i~\langle \boldsymbol{\sigma}_0 \rangle\cdot\mathbf{k} + \langle \rho_\eta \rangle~(\omega_0^j)^2~\widehat{\mathbf{u}}_0 = 0 ~. } $$ Plugging (28) into (27) gives
 * $$ \text{(30)} \qquad

\langle \boldsymbol{\sigma}_0 \rangle = \cfrac{i}{2}~ \boldsymbol{\mathsf{C}}_\text{eff} :(\widehat{\mathbf{u}}_0\otimes\mathbf{k} + \mathbf{k}\otimes\widehat{\mathbf{u}}_0)~. $$ Plugging (30) into (29) gives

-\frac{1}{2}~\left[ \boldsymbol{\mathsf{C}}_\text{eff} :(\widehat{\mathbf{u}}_0\otimes\mathbf{k} + \mathbf{k}\otimes\widehat{\mathbf{u}}_0)\right] \cdot\mathbf{k} + \langle \rho_\eta \rangle~(\omega_0^j)^2~\widehat{\mathbf{u}}_0 = 0 ~. $$ From the minor symmetry of the tensor $$\boldsymbol{\mathsf{C}}_\text{eff}$$, i.e.,

C^\text{eff}_{ijkl} = C^{\text{eff}}_{ijlk} $$ we have

\boldsymbol{\mathsf{C}}_\text{eff}:(\widehat{\mathbf{u}}_0\otimes\mathbf{k}) = \boldsymbol{\mathsf{C}}_\text{eff}:(\mathbf{k}\otimes\widehat{\mathbf{u}}_0)~. $$ Hence

\left[\boldsymbol{\mathsf{C}}_\text{eff}:(\mathbf{k}\otimes\widehat{\mathbf{u}}_0)\right]\cdot\mathbf{k} = \langle \rho_\eta \rangle~(\omega_0^j)^2~\widehat{\mathbf{u}}_0 $$ or,

{ \boldsymbol{A}_\text{eff}(\mathbf{k})~\widehat{\mathbf{u}}_0 = (\omega_0^j)^2~\widehat{\mathbf{u}}_0 } $$ where, for any vector $$\mathbf{v}$$,

{ \boldsymbol{A}_\text{eff}(\mathbf{k})~\mathbf{v} = \cfrac{1}{\langle \rho_\eta \rangle}~ \left[\boldsymbol{\mathsf{C}}_\text{eff}:(\mathbf{k}\otimes\mathbf{v})\right]\cdot\mathbf{k} ~. } $$ The quantity $$\boldsymbol{A}_\text{eff}$$ is the  effective acoustic tensor. The dispersion relations may then be calculated in a manner similar to that for electromagnetism.

Effective properties of bubbly fluid
So far we have found that, in the quasistatic limit, the Bloch solutions satisfy equations which are directly analogous to the electromagnetic or elasticity equations but with complex fields and complex effective tensors. This implies that if we have a formula for the effective tensor which is valid for real tensors, then (by analytic continuation) we can use the same formula when the tensors take complex values.

Let us examine one such situation which arises in bubbly fluids. We can assume bubbly fluids to be an assembly of coated spheres. For such a geometry, we have Hashin's relations Hashin62,Hashin62a when the material properties are real. Let us now discuss how the Hashin relations are obtained.

Hashin's relations for assemblages of coated spheres
Consider the coated sphere shown in Figure~2. The sphere has a permittivity $$\epsilon_1$$ and the coating has a permittivity $$\epsilon_2$$. Let the coated sphere be embedded in a matrix with permittivity $$\epsilon_0$$.

For certain choices of the permittivities, current will tends to flow around the coated sphere while for other choices current will be attracted toward the coated sphere. Therefore we expect that there might be a situation in which the permittivities are such that the coated sphere will have no effect on a current flowing through the matrix. In such a situation, we could continue to add spheres and completely fill space (except for a set of measure zero) without affecting the field as can be seen in Figure 3. Clearly, in that case, the effective permittivity of the assemblage of coated spheres is equal to that of the matrix.

In determining the effective properties of an assemblage of coated spheres we make the following assumptions:

 * 1) the coated spheres do not overlap the boundary of the unit cell.
 * 2) fluxes and potentials at the boundary remain unaltered due to addition of coated spheres.

The goal is to find a matrix whose properties are such that when a coated sphere is added to it, the fields in the matrix remain unaltered.

Effective Permittivity
Consider the single coated sphere shown in Figure 2. Let us assume that the sphere is centered at the origin. We look for a solution to the time-independent Maxwell's equations

\boldsymbol{\nabla} \times \mathbf{E} = 0 ~; \boldsymbol{\nabla} \cdot \mathbf{D} = 0 ~; \mathbf{D} = \epsilon~\mathbf{E} $$ with potentials

\begin{align} \varphi_1(x) & = a_1~r~\cos\theta & \qquad \text{in the core} \\ \varphi_2(x) & = \left(a_2~r + \cfrac{b_2}{r^2}\right)~\cos\theta & \qquad \text{in the coating} \\ \varphi_\text{eff}(x) & = a_\text{eff}~r~\cos\theta & \qquad \text{in the effective medium.} \end{align} $$ Then the electric field is given by

\begin{align} \mathbf{E}_1 = \boldsymbol{\nabla} \varphi_1(x) & = a_1~[\cos\theta~\mathbf{e}_r - \sin\theta~\mathbf{e}_\theta] \\ \mathbf{E}_2 = \boldsymbol{\nabla} \varphi_2(x) & = \left(a_2~r - \cfrac{2b_2}{r^3}\right)~\cos\theta~\mathbf{e}_r - \left(a_2~r + \cfrac{b_2}{r^3}\right)~\sin\theta~\mathbf{e}_\theta \\ \mathbf{E}_\text{eff} = \boldsymbol{\nabla} \varphi_\text{eff(x)} & = a_\text{eff}~ [\cos\theta~\mathbf{e}_r - \sin\theta~\mathbf{e}_\theta] ~. \end{align} $$ The potentials satisfy Laplace's equations

\nabla^2 \varphi_1 = \nabla^2 \varphi_2 = \nabla^2 \varphi_\text{eff} = 0 $$ and we only need to match the boundary conditions at the interfaces to get a solution, i.e.,
 * $$ \text{(31)} \qquad

a_1~r_c = a_2~r_c + \cfrac{b_2}{r_c^2} ~; a_\text{eff}~r_e = a_2~r_e + \cfrac{b_2}{r_e^2} ~. $$ Continuity of the tangential component of the electric field at the interface implies that
 * $$ \text{(32)} \qquad

\epsilon_1~a_1 = \epsilon_2~\left[a_2 - \cfrac{2b_2}{r_c^3}\right] ~; \epsilon_\text{eff}~a_\text{eff} = \epsilon_2~\left[a_2 - \cfrac{2b_2}{r_e^3}\right]~. $$ Combining (31) and (32) gives

\cfrac{b_2}{r_c^3} \left[1 + \cfrac{3~\epsilon_2}{\epsilon_1-\epsilon_2}\right] = \cfrac{b_2}{r_e^3} \left[1 + \cfrac{3~\epsilon_2}{\epsilon_\text{eff}-\epsilon_2}\right] ~. $$ Defining the volume fraction $$f$$ as

f_1 = 1 - f_2 = \cfrac{r_c^3}{r_e^3} $$ leads to the following expression for $$\epsilon_\text{eff}$$:

{ \epsilon_\text{eff} = \epsilon_2 + \cfrac{3~f_1~\epsilon_2(\epsilon_1-\epsilon_2)} {3~\epsilon_2 + f_2~(\epsilon_1 - \epsilon_2)} ~. } $$

Effective Bulk Modulus
In this case, consider the coated sphere shown in Figure~2 with the difference that each material now has two properties- the bulk and shear moduli. Let the bulk and shear moduli of the sphere be $$\kappa_1$$ and $$\mu_1$$, let those of the coating be $$\kappa_2$$ and $$\mu_2$$, and let the moduli of the effective medium (the matrix in the figure) be $$\kappa_\text{eff}$$ and $$\mu_\text{eff}$$, respectively. As before, the radius of the sphere is $$r_c$$ and the outer boundary of the coating has a radius $$r_e$$.

The governing equations are

\boldsymbol{\nabla} \cdot \boldsymbol{\sigma} = 0 ~; \boldsymbol{\epsilon} = \frac{1}{2}\left[\boldsymbol{\nabla}\mathbf{u} + (\boldsymbol{\nabla}\mathbf{u})^T\right] ~; \boldsymbol{\sigma} = \lambda~\text{tr}{\boldsymbol{\epsilon}}~\boldsymbol{\mathit{1}} + 2~\mu~\boldsymbol{\epsilon} ~; \kappa = \lambda + \cfrac{2}{3}~\mu ~. $$

Let the matrix be subject to a hydrostatic state of stress, i.e.,

\boldsymbol{\sigma} = -p~\boldsymbol{\mathit{1}} ~. $$ Let us look for a solution that does not perturb this field when the coated sphere is added to the matrix.

Therefore, we look for a solution with a radial displacement field

\begin{align} \mathbf{u}_1(\mathbf{x}) &= a_1~r~\mathbf{e}_r & \qquad \text{in the core} \\ \mathbf{u}_2(\mathbf{x}) &= \left(a_2~r + \cfrac{b_2}{r^2}\right)~\mathbf{e}_r & \qquad \text{in the coating} \\ \mathbf{u}_\text{eff}(\mathbf{x}) &= a_\text{eff}~r~\mathbf{e}_r & \qquad \text{in the effective medium} \end{align} $$

From the continuity of displacements at the interfaces, we have

a_1~r_c = a_2~r_c + \cfrac{b_2}{r_c^2} ~; a_\text{eff}~r_e = a_2~r_e + \cfrac{b_2}{r_e^2} ~. $$ And the continuity of radial tractions $$\boldsymbol{\sigma}\cdot\mathbf{e}_r$$ across the interface implies that

\kappa_1~a_1 = \kappa_2~a_2 - \cfrac{4~b_2~\mu_2}{3~r_c^3} ~; \kappa_\text{eff}~a_\text{eff} = \kappa_2~a_2 - \cfrac{4~b_2~\mu_2}{3~r_e^3} ~; $$ From the continuity relations and using the definition

f_1 = 1 - f_2 = \cfrac{r_c^3}{r_e^3} $$ we can show that the effective bulk modulus is given by
 * $$ \text{(33)} \qquad

{ \kappa_\text{eff} = \kappa_2 + \cfrac{f_1}{\cfrac{1}{\kappa_1 - \kappa_2} + \cfrac{f_2}{\kappa_2 + \cfrac{4}{3}~\mu_2}} ~. } $$

Bubbly Fluid
For the assemblage of coated spheres, when the bulk and shear moduli of the two phases are real, the effective bulk modulus is given by equation (33). We can then use the Bloch wave solutions to show that the same result holds when the moduli are complex.

Returning to the bubbly fluid problem, suppose that phase 1 is gas and phase 2 is water. The mixture is then a bubbly fluid. If we assume that water is incompressible, then $$\kappa_2 \rightarrow \infty$$. Hence from (33) we have
 * $$ \text{(34)} \qquad

\kappa_\text{eff} \approx \cfrac{\kappa_1 + \cfrac{4}{3}~f_2~\mu_2}{f_1} ~. $$ Also, assume that $$\kappa_1$$ (air) is independent of frequency.

Now consider a plane shear wave propagating into a bubbly fluid. Let the frequency of the wave be $$\omega$$ and let it be real. Let the wave be spatially attenuated with a complex wavevector $$\mathbf{k}$$. Then the associated strain field is given by

\boldsymbol{\epsilon}(\mathbf{x},\omega) = e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\epsilon}}(\omega) ~. $$ The stress field is given by (assuming a Newtonian fluid)

\boldsymbol{\sigma}(\mathbf{x}, \omega) = -i~\omega~\eta_\mu~e^{i~\mathbf{k}\cdot\mathbf{x}}~\widehat{\boldsymbol{\epsilon}}(\omega) = -i~\omega~\eta_\mu~\boldsymbol{\epsilon}(\mathbf{x}, \omega) $$ where $$\eta_\mu$$ is the shear viscosity of water and is assume to be independent of frequency. Therefore, we can think of water as having a complex shear modulus, i.e.,

\mu_2 = -i~\omega~\eta_\mu ~. $$ Plugging this in (34) gives

\kappa_\text{eff} \approx \cfrac{\kappa_1}{f_1} - i~\eta_\kappa^\text{eff}~\omega $$ where the effective bulk viscosity is defined as

\eta_\kappa^\text{eff} \approx \cfrac{4~f_2~\eta_\mu}{3~f_1} ~. $$ When $$f_1 \rightarrow 0$$, we get

\eta_\kappa^\text{eff} \rightarrow \cfrac{4~\eta_\mu}{3~f_1} ~. $$ Therefore, the shear viscosity of water has been converted into the bulk viscosity of the bubbly fluid. This is the reason that sound is damped strongly in bubbly fluids.


text:
Waves in composites and metamaterials/Duality relations and phase interchange identity in laminates

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.

Duality Relations in Two Dimensions
Instead of taking the direct route of the previous lecture, we can determine the effective properties of composites using duality relations.

Let us look at the quasistatic Maxwell's equations first. These equations can be written as
 * $$ \text{(1)} \qquad

\mathbf{D} = \boldsymbol{\epsilon}\cdot\mathbf{E} ~; \boldsymbol{\nabla} \cdot \mathbf{D} = 0 ~; \boldsymbol{\nabla} \times \mathbf{E} = \boldsymbol{0} $$ where $$\mathbf{D}$$ and $$\mathbf{E}$$ are periodic.

Let us define the effective permittivity of the medium ($$\boldsymbol{\epsilon}_\text{eff}$$) using
 * $$ \text{(2)} \qquad

\langle \mathbf{D} \rangle = \boldsymbol{\epsilon}_\text{eff} \cdot \langle \mathbf{E} \rangle $$ where $$\langle \mathbf{D} \rangle$$ and $$\langle \mathbf{E} \rangle$$ denote volume averages, i.e.,

\langle \mathbf{D} \rangle = \cfrac{1}{V}~\int_\Omega \mathbf{D}(\mathbf{x})~\text{d}\Omega ~; \langle \mathbf{E} \rangle = \cfrac{1}{V}~\int_\Omega \mathbf{E}(\mathbf{x})~\text{d}\Omega $$ where $$V$$ is the volume of the region $$\Omega$$.

In two dimensions, we have
 * $$ \text{(3)} \qquad

\boldsymbol{\nabla} \times \mathbf{E} = \boldsymbol{0} \qquad \implies \qquad \frac{\partial E_1}{\partial x_2} - \frac{\partial E_2}{\partial x_1} = 0 $$ since $$E_1 \equiv E_1(x_1, x_2)$$, $$E_2 \equiv E_2(x_2, x_2)$$, and $$E_3 = 0.0$$ (or constant).

Now, define
 * $$ \text{(4)} \qquad

\mathbf{D}' := \lambda~\boldsymbol{R}_\perp\cdot\mathbf{E} \qquad \text{and} \qquad \mathbf{E}' := \boldsymbol{R}_\perp \cdot \mathbf{D} $$ where $$\boldsymbol{R}_\perp$$ is the orthogonal tensor that indicates a 90$$^o$$ rotation about the $$x_3$$ axis, i.e.,

\boldsymbol{R}_\perp = \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} ~. $$ Therefore, in a rectangular Cartesian basis ($$\mathbf{e}_1, \mathbf{e}_2$$), we have

\mathbf{D}' = D'_1~\mathbf{e}_1 + D'_2~\mathbf{e}_2 = \lambda~E_2~\mathbf{e}_1 - \lambda~E_1~\mathbf{e}_2 ~; \mathbf{E}' = E'_1~\mathbf{e}_1 + E'_2~\mathbf{e}_2 = D_2~\mathbf{e}_1 - D_1~\mathbf{e}_2 ~. $$ Hence, using (3),

\boldsymbol{\nabla} \cdot \mathbf{D}' = \frac{\partial D'_1}{\partial x_1} + \frac{\partial D'_2}{\partial x_2} = \lambda\left(\frac{\partial E_2}{\partial x_1} - \frac{\partial E_1}{\partial x_2}\right) = 0 $$ and, using (1)$$_2$$,

\boldsymbol{\nabla} \times \mathbf{E}' = \left(\frac{\partial E'_2}{\partial x_1} - \frac{\partial E'_1}{\partial x_2}\right)\mathbf{e}_3 = \left(\frac{\partial D_1}{\partial x_1} + \frac{\partial D_2}{\partial x_2}\right)\mathbf{e}_3 = (\boldsymbol{\nabla} \cdot \mathbf{D})~\mathbf{e}_3 = \boldsymbol{0} ~. $$ Therefore the dual field $$\mathbf{D}'$$ and $$\mathbf{E}'$$ represent a divergence free electric field (i.e., there are no sources or sinks). Hence the field is the gradient of some potential which has zero curl.

Also, assuming that the permittivity $$\boldsymbol{\epsilon}$$ is invertible, we have from equations (4) and (1)$$_1$$

\mathbf{D}' = \lambda~\boldsymbol{R}_\perp\cdot(\boldsymbol{\epsilon}^{-1}\cdot\mathbf{D}) = \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}^{-1}\cdot (\boldsymbol{R}_\perp^T\cdot\boldsymbol{R}_\perp)\cdot\mathbf{D} = \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}^{-1}\cdot\boldsymbol{R}_\perp^T\cdot\mathbf{E}' ~. $$ Defining

\boldsymbol{\epsilon}' := \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}^{-1}\cdot\boldsymbol{R}_\perp^T $$ we then have

\mathbf{D}' = \boldsymbol{\epsilon}' \cdot \mathbf{E}' ~. $$ If we assume that $$\boldsymbol{\epsilon}$$ is symmetric, we have (in two dimensions) with respect to the basis ($$\mathbf{e}_1, \mathbf{e}_2$$)

\boldsymbol{\epsilon} \equiv \begin{bmatrix} \epsilon_{11} & \epsilon_{12} \\ \epsilon_{12} & \epsilon_{22} \end{bmatrix} \qquad \implies \qquad \boldsymbol{\epsilon}^{-1} \equiv \cfrac{1}{\det\boldsymbol{\epsilon}} \begin{bmatrix} \epsilon_{22} & -\epsilon_{12} \\ -\epsilon_{12} & \epsilon_{11} \end{bmatrix}~. $$ Then

\boldsymbol{\epsilon}' = \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}^{-1}\cdot\boldsymbol{R}_\perp^T = \cfrac{\lambda}{\det\boldsymbol{\epsilon}}~\boldsymbol{\epsilon} ~. $$ The dual system of equations is then given by
 * $$ \text{(5)} \qquad

{ \boldsymbol{\nabla} \cdot \mathbf{D}' = 0 ~; \boldsymbol{\nabla} \times \mathbf{E}' = 0 ~; \mathbf{D}' = \boldsymbol{\epsilon}'\cdot\mathbf{E}' } $$ where
 * $$ \text{(6)} \qquad

{ \mathbf{D}' = \lambda~\boldsymbol{R}_\perp\cdot\mathbf{E} ~; \mathbf{E}' = \boldsymbol{R}_\perp \cdot \mathbf{D} ~; \boldsymbol{\epsilon}' = \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}^{-1}\cdot\boldsymbol{R}_\perp^T = \cfrac{\lambda~\boldsymbol{\epsilon}}{\det\boldsymbol{\epsilon}} ~. } $$ So $$\mathbf{D}'$$ and $$\mathbf{E}'$$ solve Maxwell's equations for electricity in a  dual medium of permittivity $$\boldsymbol{\epsilon}'$$.

Effective permittivity of dual medium
The question that arises at this stage is: { what is the effective permittivity ($$\boldsymbol{\epsilon}'_\text{eff}$$) of the dual material in terms of the effective permittivity of the original material ($$\boldsymbol{\epsilon}_\text{eff}$$)?}

Taking volume averages of equations (6)$$_1$$ and (6)$$_2$$ we get

\langle \mathbf{D}' \rangle = \lambda~\boldsymbol{R}_\perp\cdot\langle \mathbf{E} \rangle ~; \langle \mathbf{E}' \rangle = \boldsymbol{R}_\perp \cdot \langle \mathbf{D} \rangle ~. $$ Recall from equation (2) that

\langle \mathbf{D} \rangle = \boldsymbol{\epsilon}_\text{eff} \cdot \langle \mathbf{E} \rangle ~. $$ Therefore,

\langle \mathbf{D}' \rangle = \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}_\text{eff}^{-1}\cdot\langle \mathbf{D} \rangle = \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}_\text{eff}^{-1} \cdot(\boldsymbol{R}^T_\perp\cdot\boldsymbol{R}_\perp)\cdot\langle \mathbf{D} \rangle = \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}_\text{eff}^{-1}\cdot\boldsymbol{R}^T_\perp\cdot\langle \mathbf{E}' \rangle~. $$ As before, defining

\boldsymbol{\epsilon}'_\text{eff} := \lambda~\boldsymbol{R}_\perp\cdot\boldsymbol{\epsilon}_\text{eff}^{-1}\cdot\boldsymbol{R}^T_\perp = \cfrac{\lambda~\boldsymbol{\epsilon}_\text{eff}}{\det(\boldsymbol{\epsilon}_\text{eff})} $$ gives

\langle \mathbf{D}' \rangle = \boldsymbol{\epsilon}'_\text{eff} \cdot \langle \mathbf{E}' \rangle ~. $$ Therefore the relation between the effective permittivity of the original and the dual material is

{ \boldsymbol{\epsilon}'_\text{eff} = \cfrac{\lambda~\boldsymbol{\epsilon}_\text{eff}}{\det(\boldsymbol{\epsilon}_\text{eff})}~. } $$

Application to a two-dimensional polycrystal
Let us consider a two-dimensional polycrystal as shown in Figure 1. The lattice vectors in each crystal are oriented in a random manner. However, for each crystal, the lattice vector can be determined by a piecewise constant rotation from a reference configuration.

Therefore, for each crystal

\boldsymbol{\epsilon}(\mathbf{x}) = \boldsymbol{R}(\mathbf{x}) \cdot \boldsymbol{\epsilon}_0 \cdot \boldsymbol{R}^T(\mathbf{x}) = \boldsymbol{R}(\mathbf{x}) \cdot \boldsymbol{\epsilon}_0 \cdot \boldsymbol{R}^{-1}(\mathbf{x}) $$ where the rotation field

\boldsymbol{R}(\mathbf{x})\cdot\boldsymbol{R}^T(\mathbf{x}) = \boldsymbol{\mathit{1}} $$ determines the local orientation of the crystal at each point and the rotations $$\boldsymbol{R}(\mathbf{x})$$ are piecewise constant in each crystal.

If we now consider a medium that is dual to the polycrystal in the sense of equation (6), then the permittivity of the dual medium is given by

\boldsymbol{\epsilon}' = \cfrac{\lambda~\boldsymbol{\epsilon}}{\det\boldsymbol{\epsilon}} ~. $$ Now

\det(\boldsymbol{\epsilon}) = \det(\boldsymbol{R})~\det(\boldsymbol{\epsilon}_0)~\det(\boldsymbol{R}^T) = \det(\boldsymbol{\epsilon}_0) $$ since

\det(\boldsymbol{R}^T) = \cfrac{1}{\det(\boldsymbol{R})} ~. $$ Therefore,

\boldsymbol{\epsilon}' = \cfrac{\lambda~\boldsymbol{\epsilon}}{\det(\boldsymbol{\epsilon}_0)} ~. $$ Let us choose the constant $$\lambda$$ such that

\lambda = \det(\boldsymbol{\epsilon}_0) ~. $$ Then

{ \boldsymbol{\epsilon}'(\mathbf{x}) = \boldsymbol{\epsilon}(\mathbf{x}) ~. } $$ Recall that

\mathbf{D} = \boldsymbol{\epsilon}\cdot\mathbf{E} \qquad \text{and} \qquad \langle \mathbf{D} \rangle = \boldsymbol{\epsilon}_\text{eff} \cdot \langle \mathbf{E} \rangle ~. $$ If

\boldsymbol{\epsilon}_\text{eff} = \langle \boldsymbol{\epsilon} \rangle \qquad \text{and} \qquad \boldsymbol{\epsilon}'_\text{eff} = \langle \boldsymbol{\epsilon}' \rangle $$ we have

\langle \boldsymbol{\epsilon}' \rangle = \langle \cfrac{\lambda~\boldsymbol{\epsilon}} {\det(\boldsymbol{\epsilon})}\rangle = \langle \boldsymbol{\epsilon} \rangle \qquad \text{if} \lambda = \det(\boldsymbol{\epsilon}) ~. $$ Hence

{ \boldsymbol{\epsilon}'_\text{eff} = \boldsymbol{\epsilon}_\text{eff} ~. } $$ But we also have

\boldsymbol{\epsilon}'_\text{eff} = \cfrac{\lambda}{\det(\boldsymbol{\epsilon}_\text{eff})}~\boldsymbol{\epsilon}_\text{eff}. $$ Therefore,

{ \lambda = \det(\boldsymbol{\epsilon}_\text{eff}) = \det(\boldsymbol{\epsilon}_0)~. } $$ In particular, if the polycrystal is isotropic, i.e., $$\boldsymbol{\epsilon}_\text{eff} = \epsilon_\text{eff}~\boldsymbol{\mathit{1}}$$, then we have ({\Red Show this!})

{ \epsilon_\text{eff} = \sqrt{\det(\boldsymbol{\epsilon}_0)} } $$ In this case only the root with the positive imaginary part is the correct solution unless $$\det{\boldsymbol{\epsilon}_0}$$ is real in which case only the positive root is correct.

Application to a 2-D composite of two isotropic phases
Consider the composite of two isotropic phases shown in Figure 2. Define an indicator function

\chi(\mathbf{x}) = \begin{cases} 1 & \text{in phase}~1 \\ 0 & \text{otherwise} \end{cases} $$

Then, since the phases are isotropic, we can write

\boldsymbol{\epsilon}(\mathbf{x}) = \left[\epsilon_1~\chi(\mathbf{x}) + \epsilon_2~\{1 - \chi(\mathbf{x})\} \right]~\boldsymbol{\mathit{1}} \qquad \implies \qquad \det(\boldsymbol{\epsilon}(\mathbf{x})) = \left[\epsilon_1~\chi(\mathbf{x}) + \epsilon_2~\{1 - \chi(\mathbf{x})\} \right]^2 ~. $$ The dual material is defined as one having a permittivity given by

\boldsymbol{\epsilon}'(\mathbf{x}) = \cfrac{\lambda~\boldsymbol{\epsilon}(\mathbf{x})}{\det(\boldsymbol{\epsilon}(\mathbf{x}))} = \cfrac{\lambda}{\epsilon_1~\chi(\mathbf{x}) + \epsilon_2~\{1 - \chi(\mathbf{x})\}}~ \boldsymbol{\mathit{1}} ~. $$ We can write the above as

\boldsymbol{\epsilon}'(\mathbf{x}) = \begin{cases} \cfrac{\lambda}{\epsilon_1}~\boldsymbol{\mathit{1}} & \text{in phase} ~1 \\ \cfrac{\lambda}{\epsilon_2}~\boldsymbol{\mathit{1}} & \text{otherwise} \end{cases} $$ or alternatively,

\boldsymbol{\epsilon}'(\mathbf{x}) = \left[\cfrac{\lambda}{\epsilon_1}~\chi(\mathbf{x}) + \cfrac{\lambda}{\epsilon_2}~\{1 - \chi(\mathbf{x})\}\right]~\boldsymbol{\mathit{1}} ~. $$ If we choose $$\lambda = \epsilon_1~\epsilon_2$$, we get

\boldsymbol{\epsilon}'(\mathbf{x}) = \left[\epsilon_2~\chi(\mathbf{x}) + \epsilon_1~\{1 - \chi(\mathbf{x})\}\right]~\boldsymbol{\mathit{1}} ~. $$  Note that the phases are interchanged in the dual material!

Now, recall that the effective permittivities of the original and the dual material are related by

\boldsymbol{\epsilon}'_\text{eff} = \cfrac{\lambda~\boldsymbol{\epsilon}_\text{eff}}{\det(\boldsymbol{\epsilon}_\text{eff})} = \cfrac{\epsilon_1~\epsilon_2~\boldsymbol{\epsilon}_\text{eff}}{\det(\boldsymbol{\epsilon}_\text{eff})} ~. $$

We can use this relation to find the effective permittivities of materials that are invariant with respect to interchange of phases. Examples of such materials are checkerboard material and random polycrystals where each crystal has an equal probability of being of phase 1 or phase 2.

For such a phase interchange invariant material, the effective permittivity of the original material is equal to that of the dual material, i.e.,

\boldsymbol{\epsilon}'_\text{eff} = \boldsymbol{\epsilon}_\text{eff} \qquad \implies \qquad \det(\boldsymbol{\epsilon}_\text{eff}) = \epsilon_1~\epsilon_2 ~. $$ If the composite material is isotropic, i.e., $$\boldsymbol{\epsilon}_\text{eff} = \epsilon_\text{eff}~\boldsymbol{\mathit{1}}$$, then

\det(\boldsymbol{\epsilon}_\text{eff}) = \epsilon_\text{eff}^2 ~. $$ Hence,

{ \epsilon_\text{eff} = \sqrt{\epsilon_1~\epsilon_2} ~. } $$ This is an useful result that can be used to test numerical codes.

Paradox
If $$\epsilon_1 = 1$$ and $$\epsilon_2 = -1$$ then both materials are lossless. \footnote{\Red Need to add section here showing that energy dissipation is proportional to $$\text{Im}(\boldsymbol{\epsilon})$$.} But $$\epsilon_\text{eff} = \sqrt{-1} = i$$. So the composite dissipates energy into heat. But where?

To see this we should take $$\epsilon_2 = -1 + i\delta$$ and look at the limit where $$\delta \rightarrow 0$$. In this limit, the fields lose their square integrability at the corners (in a checkerboard). So an enormous amount of heat per unit volume is dissipated in the vicinity of each corner.

Extensions to 2-D elasticity
Duality and phase interchange relations for elasticity were first derived by Berdichevski~Berdi83. In that work, an exact formula for the shear modulus of a checkerboard material (with two incompressible phases) was derived. Further extensions and details of can be found in Sections 3.5, 3.6, 3.7 and 4.7 in Milton02.

We can apply duality transformations to incompressible media or media where the bulk modulus $$\kappa$$ is equal in both phases and the shear moduli of the two phases are $$\mu_1$$ and $$\mu_2$$. For example, if we have a phase interchange invariant composite that is isotropic and two-dimensional (such as a checkerboard or a cell material), then the effective elastic moduli are given by

\begin{align} \kappa_\text{eff} & = \kappa \\ \mu_\text{eff} & = \cfrac{\kappa}{-1+ \sqrt{\left(1 + \cfrac{\kappa}{\mu_1}\right) \left(1 + \cfrac{\kappa}{\mu_2}\right)}}~. \end{align} $$ In the case where $$\mu_1 = \mu_2$$ we have $$\mu_\text{eff} = \mu$$.

The Effective Tensors of Laminate Materials
In this section we will discuss the method of Backus~Backus62. Similar approaches have also been used by Postma~Postma55 and Tartar~Tartar76.

Consider a material laminated in the $$x_1$$ direction as shown in Figure 3.

To find the relation between $$\langle \mathbf{D} \rangle$$ and $$\langle \mathbf{E} \rangle$$ we cannot average the constitutive relation

\mathbf{D} = \boldsymbol{\epsilon}\cdot\mathbf{E} \qquad \implies \qquad \langle \mathbf{D} \rangle = \langle \boldsymbol{\epsilon}\cdot\mathbf{E} \rangle $$ because

\langle \boldsymbol{\epsilon}\cdot\mathbf{E} \rangle \ne \langle \boldsymbol{\epsilon} \rangle\cdot\langle \mathbf{E} \rangle $$ unless $$\boldsymbol{\epsilon}$$ is constant or $$\mathbf{E}$$ is constant.

However, there are fields which are constant in certain directions and those can be used to simplify things. Since the tangential components (parallel to the layers) of the electric field ($$\mathbf{E}$$) are piecewise constant and continuous across the interfaces between the layers, these tangential components must be constant, i.e., $$E_2$$ and $$E_3$$ are constant in the laminate. Similarly, the continuity of the normal electric displacement field ($$\mathbf{D}$$) across the interfaces and the fact that this field is constant in each layer implies that the component $$D_1$$ is constant in the laminate.

Let us rewrite the constitutive relation in matrix form (with respect to the rectangular Cartesian basis ($$\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$$)) so that constant fields appear on the right hand side, i.e.,
 * $$ \text{(7)} \qquad

\begin{bmatrix} D_1 \\ \mathbf{D}_t \end{bmatrix} = \begin{bmatrix} \epsilon_{11} & \boldsymbol{\epsilon}_{1t} \\ \boldsymbol{\epsilon}_{t1}   & \boldsymbol{\epsilon}_{tt} \end{bmatrix} \begin{bmatrix} E_1 \\ \mathbf{E}_t \end{bmatrix} $$ where
 * $$ \text{(8)} \qquad

\mathbf{D}_t = \begin{bmatrix} D_2 \\ D_3 \end{bmatrix} ~; \mathbf{E}_t = \begin{bmatrix} E_2 \\ E_3 \end{bmatrix} ~; \boldsymbol{\epsilon}_{1t} = \begin{bmatrix} \epsilon_{12} & \epsilon_{13}\end{bmatrix} ~; \boldsymbol{\epsilon}_{t1} = \begin{bmatrix} \epsilon_{21} \\ \epsilon_{31}\end{bmatrix}~; \boldsymbol{\epsilon}_{tt} = \begin{bmatrix} \epsilon_{22} & \epsilon_{23}\\ \epsilon_{32} & \epsilon_{33} \end{bmatrix}~. $$ Note that the constant fields are $$D_1$$ and $$\mathbf{E}_t$$. We want to rewrite equation (7) so that these constant fields appear on the right hand side.

From the first row of (7) we get

D_1 = \epsilon_{11}~E_1 + \boldsymbol{\epsilon}_{1t}\cdot\mathbf{E}_t $$ or,
 * $$ \text{(9)} \qquad

E_1 = \cfrac{1}{\epsilon_{11}}~D_1 - \cfrac{\boldsymbol{\epsilon}_{1t}}{\epsilon_{11}} \cdot\mathbf{E}_t ~. $$ From the second row of (7) we get
 * $$ \text{(10)} \qquad

\mathbf{D}_t = \boldsymbol{\epsilon}_{t1}~E_1 + \boldsymbol{\epsilon}_{tt}\cdot\mathbf{E}_t ~. $$ Substitution of (9) into (10) gives
 * $$ \text{(11)} \qquad

\mathbf{D}_t = \cfrac{\boldsymbol{\epsilon}_{t1}}{\epsilon_{11}}~D_1 - \cfrac{\boldsymbol{\epsilon}_{t1}\cdot\boldsymbol{\epsilon}_{1t}}{\epsilon_{11}}\cdot\mathbf{E}_t + \boldsymbol{\epsilon}_{tt}\cdot\mathbf{E}_t ~. = \cfrac{\boldsymbol{\epsilon}_{t1}}{\epsilon_{11}}~D_1 + \left(\boldsymbol{\epsilon}_{tt}          - \cfrac{\boldsymbol{\epsilon}_{t1}\cdot\boldsymbol{\epsilon}_{1t}}{\epsilon_{11}}\right) \cdot\mathbf{E}_t ~. $$ Collecting (9) and (11) gives

\begin{bmatrix} -E_1 \\ \mathbf{D}_t \end{bmatrix} = \cfrac{1}{\epsilon_{11}}~\begin{bmatrix} -1 & \boldsymbol{\epsilon}_{1t} \\ \boldsymbol{\epsilon}_{t1} & \epsilon_{11}~\boldsymbol{\epsilon}_{tt} - \boldsymbol{\epsilon}_{t1}\cdot\boldsymbol{\epsilon}_{1t} \end{bmatrix} \begin{bmatrix} -D_1 \\ \mathbf{E}_t \end{bmatrix} $$ where the negative signs on $$E_1$$ and $$D_1$$ are used to make sure that the signs of the off diagonal terms are identical.

Define

{ \boldsymbol{L}(\mathbf{x}) := \cfrac{1}{\epsilon_{11}}~\begin{bmatrix} -1 & \boldsymbol{\epsilon}_{1t} \\ \boldsymbol{\epsilon}_{t1} & \epsilon_{11}~\boldsymbol{\epsilon}_{tt} - \boldsymbol{\epsilon}_{t1}\cdot\boldsymbol{\epsilon}_{1t} \end{bmatrix} ~. } $$ Then we have,
 * $$ \text{(12)} \qquad

{ \begin{bmatrix} -E_1 \\ \mathbf{D}_t \end{bmatrix} = \boldsymbol{L}(\mathbf{x})\cdot \begin{bmatrix} -D_1 \\ \mathbf{E}_t \end{bmatrix} ~. } $$ Since the vector on the right hand side is constant, an volume average of (12) gives
 * $$ \text{(13)} \qquad

\begin{bmatrix} -\langle E_1 \rangle \\ \langle \mathbf{D}_t \rangle \end{bmatrix} = \langle \boldsymbol{L}(\mathbf{x}) \rangle\cdot \begin{bmatrix} -D_1 \\ \mathbf{E}_t \end{bmatrix} ~. $$

Let us define the effective permittivity of the laminate $$\boldsymbol{\epsilon}_\text{eff}$$ via

\langle \mathbf{D} \rangle = \boldsymbol{\epsilon}_\text{eff} \cdot \langle \mathbf{E} \rangle ~. $$ Since the tangential components of $$\mathbf{E}$$ are constant in the laminate, the average values $$\langle E_2 \rangle$$ and $$\langle {E_3}\rangle$$ must also be constant. Similarly, the average value $$\langle D_1 \rangle$$ must be constant. Therefore we can use the same arguments as we used before to write the effective constitutive relation in the form
 * $$ \text{(14)} \qquad

{ \begin{bmatrix} -\langle E_1 \rangle \\ \langle \mathbf{D}_t \rangle \end{bmatrix} = \boldsymbol{L}_\text{eff} \cdot \begin{bmatrix} -\langle D_1 \rangle \\ \langle \mathbf{E}_t \rangle\end{bmatrix} = \boldsymbol{L}_\text{eff} \cdot \begin{bmatrix} -D_1 \\ \mathbf{E}_t\end{bmatrix} } $$ where

{ \boldsymbol{L}_\text{eff} := \cfrac{1}{\epsilon^\text{eff}_{11}}~\begin{bmatrix} -1 & \boldsymbol{\epsilon}^\text{eff}_{1t} \\ \boldsymbol{\epsilon}^\text{eff}_{t1} & \epsilon_{11}~\boldsymbol{\epsilon}^\text{eff}_{tt} - \boldsymbol{\epsilon}^\text{eff}_{t1}\cdot\boldsymbol{\epsilon}_{1t} \end{bmatrix} } $$ and $$\boldsymbol{\epsilon}_\text{eff}$$ has been decomposed in exactly the same manner as $$\boldsymbol{\epsilon}$$ (see equation (8).

If we compare equations (13) and (14) we get

{ \boldsymbol{L}_\text{eff} = \langle \boldsymbol{L}(\mathbf{x}) \rangle ~. } $$ Thus we have a formula for determining the effective permittivity of the laminate.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Waves in composites and metamaterials/Backus formula for laminates and rank-1 laminates

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
Recall the material laminated in the $$x_1$$ direction as shown in Figure 1.

We showed that we could write

\begin{bmatrix} -E_1 \\ \mathbf{D}_t \end{bmatrix} = \boldsymbol{L}(x_1)\cdot \begin{bmatrix} -D_1 \\ \mathbf{E}_t \end{bmatrix} $$ or,

\begin{bmatrix} -\langle E_1 \rangle \\ \langle \mathbf{D}_t \rangle \end{bmatrix} = \langle \boldsymbol{L}(x_1) \rangle\cdot \begin{bmatrix} -D_1 \\ \mathbf{E}_t \end{bmatrix} $$ where

\boldsymbol{L}(x_1) := \cfrac{1}{\epsilon_{11}}~\begin{bmatrix} -1 & \boldsymbol{\epsilon}_{1t} \\ \boldsymbol{\epsilon}_{t1} & \epsilon_{11}~\boldsymbol{\epsilon}_{tt} - \boldsymbol{\epsilon}_{t1}\cdot\boldsymbol{\epsilon}_{1t} \end{bmatrix} ~. $$ We also showed that

\begin{bmatrix} -\langle E_1 \rangle \\ \langle \mathbf{D}_t \rangle \end{bmatrix} = \boldsymbol{L}_\text{eff} \cdot \begin{bmatrix} -D_1 \\ \mathbf{E}_t\end{bmatrix} $$ where

\boldsymbol{L}_\text{eff} := \cfrac{1}{\epsilon^\text{eff}_{11}}~\begin{bmatrix} -1 & \boldsymbol{\epsilon}^\text{eff}_{1t} \\ \boldsymbol{\epsilon}^\text{eff}_{t1} & \epsilon_{11}~\boldsymbol{\epsilon}^\text{eff}_{tt} - \boldsymbol{\epsilon}^\text{eff}_{t1}\cdot\boldsymbol{\epsilon}_{1t} \end{bmatrix} ~. $$ Therefore, we got the relation

\boldsymbol{L}_\text{eff} = \langle \boldsymbol{L}(x_1) \rangle ~. $$ This provides the relations

{ \begin{align} \epsilon_{11}^\text{eff} & = \langle \cfrac{1}{\epsilon_{11}} \rangle^{-1} \\ \epsilon_{1j}^\text{eff} & = \langle \cfrac{1}{\epsilon_{11}} \rangle^{-1} \langle \cfrac{\epsilon_{1j}}{\epsilon_{11}} \rangle \\ \epsilon_{ij}^\text{eff} & = \langle \epsilon_{ij} - \cfrac{\epsilon_{i1}~\epsilon_{1j}}{\epsilon_{11}}\rangle + \langle \cfrac{1}{\epsilon_{11}} \rangle^{-1} \langle \cfrac{\epsilon_{i1}}{\epsilon_{11}} \rangle \langle \cfrac{\epsilon_{1j}}{\epsilon_{11}} \rangle, i,j \ne 1 \end{align} } $$ where

\langle f(x_1) \rangle = \cfrac{1}{V}~\int_0^l f(x_1)~\text{d}x_1 ~. $$ When the off diagonal elements vanish, we get

\begin{align} \epsilon_{11}^\text{eff} & = \langle \cfrac{1}{\epsilon_{11}}\rangle^{-1} & & \qquad \text{(Harmonic average)} \\ \epsilon_{jj}^\text{eff} & = \langle \epsilon_{jj \rangle} & & \qquad \text{(Arithmetic average)} \\ \epsilon_{ab}^\text{eff} & = 0 & & \qquad \text{if}~ a \ne b ~. \end{align} $$ The harmonic average corresponds to a situation in which each layer may be thought of as a capacitor in series while the arithmetic average corresponds to a situation where the capacitors are in parallel.

Effective Elastic Properties of a Layered Medium
In this section we use the approach of arranging constant fields to find the effective elastic properties of a layered medium in which each layer is anisotropic. The approach used is that of Backus (Backus62).

Consider the layered medium shown in Figure 2. In this case the displacement field is continuous across the interfaces between the layers.

The strain field ($$\boldsymbol{\varepsilon}$$) is given by

\boldsymbol{\varepsilon}(\mathbf{x}) = \frac{1}{2} \left[\boldsymbol{\nabla}\mathbf{u} + (\boldsymbol{\nabla}\mathbf{u})^T\right] ~. $$ In a Cartesian basis ($$\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$$) with coordinates ($$x_1, x_2, x_3$$) we can write

\varepsilon_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) ~. $$ Therefore the strain components $$\varepsilon_{22}, \varepsilon_{33}, \varepsilon_{23}$$ are also continuous across the interfaces. Moreover, these components of strain are also constant in each layer. Since a piecewise constant field that is also continuous must be constant, the strain components $$\varepsilon_{22}, \varepsilon_{33}, \varepsilon_{23}$$ must be constant throughout the laminate.

The tractions (normal components of the stress) at each interface are given by

\mathbf{t} = \boldsymbol{\sigma}\cdot\mathbf{n} \qquad \implies \qquad t_i = \sigma_{ij}~n_j $$ where $$\boldsymbol{\sigma}$$ is the stress and $$\mathbf{n}$$ is the normal at the interface. Now the tractions must be continuous at the interfaces. Since the normal components of the stress are piecewise constant in each layer, this implies that the normal components of the stress must also be constant throughout the laminate.

We have chosen $$\mathbf{n}$$ such that $$\mathbf{n} = (1, 0, 0)$$. Therefore the stress components $$\sigma_{11}, \sigma_{12}, \sigma_{13}$$ must be constant.

Recall that the constitutive relation for an anisotropic elastic material is given by

\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}}:\boldsymbol{\varepsilon} ~. $$ Following the approach that we used for the permittivity, we now write the constitutive relation in the form

\begin{bmatrix} \boldsymbol{\sigma}_n \\ \boldsymbol{\sigma}_t \end{bmatrix} = \begin{bmatrix} \boldsymbol{C}_{nn} & \boldsymbol{C}_{nt} \\ \boldsymbol{C}_{tn} & \boldsymbol{C}_{tt} \end{bmatrix} \begin{bmatrix} \boldsymbol{\varepsilon}_n \\ \boldsymbol{\varepsilon}_t \end{bmatrix} $$ where

\boldsymbol{\sigma}_n = \begin{bmatrix} \sigma_{11} \\ \sqrt{2}~\sigma_{12} \\ \sqrt{2}~\sigma_{13} \end{bmatrix} ~; \boldsymbol{\sigma}_t = \begin{bmatrix} \sigma_{22} \\ \sigma_{33} \\ \sqrt{2}~\sigma_{23} \end{bmatrix} ~; \boldsymbol{\varepsilon}_n = \begin{bmatrix} \varepsilon_{11} \\ \sqrt{2}~\varepsilon_{12} \\ \sqrt{2}~\varepsilon_{13} \end{bmatrix} ~; \boldsymbol{\varepsilon}_t = \begin{bmatrix} \varepsilon_{22} \\ \varepsilon_{33} \\ \sqrt{2}~\varepsilon_{23} \end{bmatrix} $$ and

\boldsymbol{C}_{nn} = \begin{bmatrix} C_{1111}         & \sqrt{2}~C_{1112} & \sqrt{2}~C_{1113}\\ \sqrt{2}~C_{1211} & 2~C_{1212}       & 2~C_{1213}\\ \sqrt{2}~C_{1311} & 2~C_{1312}       & 2~C_{1313} \end{bmatrix} ~; \boldsymbol{C}_{nt} = \begin{bmatrix} C_{1122}         & C_{1133}          & \sqrt{2}~C_{1123}\\ \sqrt{2}~C_{1222} & \sqrt{2}~C_{1233} & 2~C_{1223}\\ \sqrt{2}~C_{1322} & \sqrt{2}~C_{1333} & 2~C_{1323} \end{bmatrix} ~; $$

\boldsymbol{C}_{tn} = \begin{bmatrix} C_{2211}         & \sqrt{2}~C_{2212} & \sqrt{2}~C_{2213}\\ C_{3311}         & \sqrt{2}~C_{3312} & \sqrt{2}~C_{3313}\\ \sqrt{2}~C_{2311} & 2~C_{2312}       & 2~C_{2313} \end{bmatrix} ~; \boldsymbol{C}_{tt} = \begin{bmatrix} C_{2222}         & C_{2233}          & \sqrt{2}~C_{2223}\\ C_{3322}         & C_{3333}          & \sqrt{2}~C_{3323}\\ \sqrt{2}~C_{2322} & \sqrt{2}~C_{2333} & 2~C_{2323} \end{bmatrix} ~. $$ From the major symmetry of $$\boldsymbol{\mathsf{C}}$$, we see that $$\boldsymbol{C}_{nt} = \boldsymbol{C}_{tn}^T$$. Also, $$\boldsymbol{C}_{nn}$$ and $$\boldsymbol{C}_{tt}$$ are symmetric.

Writing the first row out, we get
 * $$ \text{(1)} \qquad

\boldsymbol{\sigma}_n = \boldsymbol{C}_{nn}\cdot\boldsymbol{\varepsilon}_n + \boldsymbol{C}_{nt}\cdot\boldsymbol{\varepsilon}_t \qquad \implies \qquad \boldsymbol{\varepsilon}_n = \boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{\sigma}_n - \left[\boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt}\right]\cdot\boldsymbol{\varepsilon}_t~. $$ From the second row we get

\boldsymbol{\sigma}_t = \boldsymbol{C}_{nt}^T\cdot\boldsymbol{\varepsilon}_n + \boldsymbol{C}_{tt}\cdot\boldsymbol{\varepsilon}_t ~. $$ Substituting the expression for $$\boldsymbol{\varepsilon}_n$$ from the first row, we get
 * $$ \text{(2)} \qquad

\boldsymbol{\sigma}_t = \boldsymbol{C}_{nt}^T\cdot\left[\boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{\sigma}_n - \boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt}\cdot\boldsymbol{\varepsilon}_t\right] + \boldsymbol{C}_{tt}\cdot\boldsymbol{\varepsilon}_t = \left[\boldsymbol{C}_{nt}^T\cdot\boldsymbol{C}_{nn}^{-1}\right]\cdot\boldsymbol{\sigma}_n + \left[\boldsymbol{C}_{tt} - \boldsymbol{C}_{nt}^T\cdot\boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt}\right] \cdot\boldsymbol{\varepsilon}_t ~. $$ Collecting (1) and (2) we get

\begin{bmatrix} -\boldsymbol{\varepsilon}_n \\ \\\boldsymbol{\sigma}_t \end{bmatrix} = \begin{bmatrix} -\boldsymbol{C}_{nn}^{-1} & \boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt} \\ \\   \boldsymbol{C}_{nt}^T\cdot\boldsymbol{C}_{nn}^{-1} & \boldsymbol{C}_{tt} - \boldsymbol{C}_{nt}^T\cdot\boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt} \end{bmatrix} \begin{bmatrix} -\boldsymbol{\sigma}_n \\ \\ \boldsymbol{\varepsilon}_t \end{bmatrix} ~. $$ Taking a volume average gives
 * $$ \text{(3)} \qquad

\begin{bmatrix} -\langle \boldsymbol{\varepsilon}_n \rangle \\ \\\langle \boldsymbol{\sigma}_t \rangle \end{bmatrix} = \begin{bmatrix} -\langle \boldsymbol{C}_{nn}^{-1} \rangle & \langle \boldsymbol{C}_{nn} ^{-1}\cdot\boldsymbol{C}_{nt}\rangle \\ \\   \langle \boldsymbol{C}_{nt}^T\cdot\boldsymbol{C}_{nn}^{-1} \rangle & \langle \boldsymbol{C}_{tt} - \boldsymbol{C}_{nt}^T\cdot\boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt}\rangle \end{bmatrix} \begin{bmatrix} -\boldsymbol{\sigma}_n \\ \\ \boldsymbol{\varepsilon}_t \end{bmatrix} ~. $$ If the effective stiffness of the material is given by

\langle \boldsymbol{\sigma} \rangle = \boldsymbol{\mathsf{C}}_\text{eff}:\langle \boldsymbol{\varepsilon} \rangle $$ we can also show that
 * $$ \text{(4)} \qquad

\begin{bmatrix} -\langle \boldsymbol{\varepsilon}_n \rangle \\ \\\langle \boldsymbol{\sigma}_t \rangle \end{bmatrix} = \begin{bmatrix} -(\boldsymbol{C}^\text{eff}_{nn})^{-1} & (\boldsymbol{C}^\text{eff}_{nn})^{-1}\cdot\boldsymbol{C}^\text{eff}_{nt} \\ \\   (\boldsymbol{C}^\text{eff}_{nt})^T\cdot(\boldsymbol{C}^\text{eff}_{nn})^{-1} & \boldsymbol{C}^\text{eff}_{tt} - (\boldsymbol{C}^\text{eff}_{nt})^T\cdot(\boldsymbol{C}^\text{eff}_{nn})^{-1} \cdot\boldsymbol{C}^\text{eff}_{nt} \end{bmatrix} \begin{bmatrix} -\boldsymbol{\sigma}_n \\ \\ \boldsymbol{\varepsilon}_t \end{bmatrix} ~. $$ Comparing (3) and (4) we can show that

{ \begin{align} \boldsymbol{C}_{nn}^\text{eff} & = \left\langle \boldsymbol{C}_{nn}^{-1} \right\rangle^{-1} \\ \boldsymbol{C}_{nt}^\text{eff} & = \left\langle \boldsymbol{C}_{nn}^{-1} \right\rangle^{-1}\cdot \left\langle \boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt} \right\rangle \\ \boldsymbol{C}_{tt}^\text{eff} & = \left\langle \boldsymbol{C}_{tt} - \boldsymbol{C}_{tn}\cdot\boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt} \right\rangle + \left\langle \boldsymbol{C}_{tn}\cdot\boldsymbol{C}_{nn}^{-1} \right\rangle\cdot \left\langle \boldsymbol{C}_{nn}^{-1} \right\rangle^{-1}\cdot \left\langle \boldsymbol{C}_{nn}^{-1}\cdot\boldsymbol{C}_{nt} \right\rangle~. \end{align} } $$

Isotropic layers
If the material in each layer is isotropic, then the constitutive relation is

\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}}:\boldsymbol{\epsilon} = \lambda(x_1)~\text{tr}(\boldsymbol{\varepsilon})~\boldsymbol{\mathit{1}} + 2~\mu(x_1)~\boldsymbol{\varepsilon} $$ where $$\lambda$$ is the Lame modulus and $$\mu$$ is the shear modulus. In that case the effective properties of the laminate are

{ \begin{align} C_{1111}^\text{eff} & = \langle \cfrac{1}{\lambda+2~\mu}\rangle^{-1} ~; C_{1122}^\text{eff} = C_{1133}^\text{eff} = \langle \cfrac{\lambda}{\lambda+2~\mu} \rangle~\langle \cfrac{1}{\lambda+2~\mu}\rangle^{-1}\\ C_{1212}^\text{eff} & = C_{1313}^\text{eff} = \langle \cfrac{1}{\mu} \rangle^{-1} ~; C_{2323}^\text{eff} = \langle \mu \rangle ~; C_{1112}^\text{eff} = C_{1113}^\text{eff} = C_{1123}^\text{eff} = 0 \\ C_{2222}^\text{eff} & = C_{3333}^\text{eff} = \langle \cfrac{4~\mu~\lambda~(\lambda+\mu)}{(\lambda+2~\mu)^2} \rangle + \langle \cfrac{\lambda}{\lambda+2~\mu} \rangle~\langle \cfrac{1}{\lambda+2~\mu}\rangle^{-1} \\ C_{2233}^\text{eff} & = \langle \cfrac{2~\mu~\lambda}{\lambda+2~\mu} \rangle + \langle \cfrac{\lambda}{\lambda+2~\mu}\rangle^2~\langle \cfrac{1} {\lambda+2~\mu}\rangle^{-1} \end{align} } $$

Laminates with Arbitrary Direction of Lamination n
So far we have dealt with laminates with a single direction of lamination that was oriented with the $$x_1$$ axis. In this section we generalize our approach to deal with laminates with a normal $$\mathbf{n}$$ which is not parallel to the $$x_1$$ axis.

Recall that the normal component of $$\mathbf{D}$$, i.e., $$\mathbf{D}\cdot\mathbf{n}$$, is constant and the tangential components of $$\mathbf{E}$$ are constant throughout the entire laminate (if there is only one direction of lamination).

Let us introduce the second-order tensor basis

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) = \mathbf{n}\otimes\mathbf{n} ~; \boldsymbol{\mathit{\Gamma}}_2(\mathbf{n}) = \boldsymbol{\mathit{1}} - \mathbf{n}\otimes\mathbf{n} = \boldsymbol{\mathit{1}} - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) ~. $$

These are useful because
 * $$ \text{(5)} \qquad

\begin{align} \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{D}(\mathbf{x}) & = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{D} \rangle & & \quad \text{normal component} \\ \boldsymbol{\mathit{\Gamma}}_2(\mathbf{n})\cdot\mathbf{E}(\mathbf{x}) & = \boldsymbol{\mathit{\Gamma}}_2(\mathbf{n})\cdot\langle \mathbf{E} \rangle & & \quad \text{tangential components.} \end{align} $$ Therefore,
 * $$ \text{(6)} \qquad

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{E}(\mathbf{x}) = \mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle + \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{E} \rangle ~. $$ Let us now introduce a polarization field
 * $$ \text{(7)} \qquad

\mathbf{P}(\mathbf{x}) = [\boldsymbol{\epsilon}(\mathbf{x}) - \epsilon_0~\boldsymbol{\mathit{1}}]\cdot\mathbf{E}(\mathbf{x}) = \mathbf{D}(\mathbf{x}) - \epsilon_0~\mathbf{E}(\mathbf{x}) $$ where $$\epsilon_0$$ is an arbitrary constant.

The volume averaged polarization field is given by
 * $$ \text{(8)} \qquad

\langle \mathbf{P} \rangle = \boldsymbol{\epsilon}_\text{eff}\cdot\langle \mathbf{E} \rangle - \epsilon_0~\langle \mathbf{E} \rangle = (\boldsymbol{\epsilon}_\text{eff} - \epsilon_0~\boldsymbol{\mathit{1}})\cdot\langle \mathbf{E} \rangle ~. $$ Define

\begin{align} \boldsymbol{S}(\mathbf{x}) & := \epsilon_0~[\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}(\mathbf{x})]^{-1} \\ \boldsymbol{S}_\text{eff} & := \epsilon_0~[\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]^{-1}~. \end{align} $$ Then,
 * $$ \text{(9)} \qquad

\begin{align} \boldsymbol{S}(\mathbf{x})\cdot\mathbf{P}(\mathbf{x}) & = -\left(\epsilon_0~[\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}(\mathbf{x})]^{-1}\right) \cdot\left( [\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}(\mathbf{x})]\cdot\mathbf{E}(\mathbf{x})\right) = -\epsilon_0~\mathbf{E}(\mathbf{x}) \\ \boldsymbol{S}_\text{eff}\cdot\langle \mathbf{P} \rangle & = -\left(\epsilon_0~[\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]^{-1}\right) \cdot\left([\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]\cdot\langle \mathbf{E} \rangle\right) = -\epsilon_0~\langle \mathbf{E} \rangle ~. \end{align} $$ Applying the projection $$\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})$$ to (7), we get

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{P}(\mathbf{x}) = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{D}(\mathbf{x}) - \epsilon_0~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{E}(\mathbf{x}) ~. $$ Using equations (5)$$_1$$ and (6), we have

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{P}(\mathbf{x}) = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{D} \rangle - \epsilon_0~\mathbf{E}(\mathbf{x}) + \epsilon_0~\langle \mathbf{E} \rangle - \epsilon_0~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{E} \rangle ~. $$ From the definitions (9) we can then write

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{P}(\mathbf{x}) = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{D} \rangle + \boldsymbol{S}(\mathbf{x})\cdot\mathbf{P}(\mathbf{x}) - \boldsymbol{S}_\text{eff}\cdot\langle \mathbf{P} \rangle - \epsilon_0~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{E} \rangle ~. $$ Define
 * $$ \text{(10)} \qquad

\mathbf{V}(\mathbf{n},\mathbf{x}) := \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{D} \rangle - \boldsymbol{S}_\text{eff}\cdot\langle \mathbf{P} \rangle - \epsilon_0~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{E} \rangle ~. $$ Then we have

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{P}(\mathbf{x}) = \boldsymbol{S}(\mathbf{x})\cdot\mathbf{P}(\mathbf{x}) + \mathbf{V}(\mathbf{n}, \mathbf{x}) $$ or
 * $$ \text{(11)} \qquad

\mathbf{V}(\mathbf{n}, \mathbf{x}) = -\left[\boldsymbol{S}(\mathbf{x}) - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]\cdot\mathbf{P}(\mathbf{x})~. $$ Also, form equations (10) and (8) we have

\mathbf{V}(\mathbf{n},\mathbf{x}) = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\left[\langle \mathbf{D} \rangle - \epsilon_0~\cdot\langle \mathbf{E} \rangle\right] - \boldsymbol{S}_\text{eff}\cdot\langle \mathbf{P} \rangle = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{P} \rangle - \boldsymbol{S}_\text{eff}\cdot\langle \mathbf{P} \rangle $$ or,
 * $$ \text{(12)} \qquad

\mathbf{V}(\mathbf{n}, \mathbf{x}) = -\left[\boldsymbol{S}_\text{eff} - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]\cdot\langle \mathbf{P} \rangle ~. $$ Inverting (11) and (12) we have
 * $$ \text{(13)} \qquad

\begin{align} \mathbf{P}(\mathbf{x}) & = - \left[\boldsymbol{S}(\mathbf{x}) - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1}\cdot \mathbf{V}(\mathbf{n}, \mathbf{x}) \\ \langle \mathbf{P} \rangle & = -\left[\boldsymbol{S}_\text{eff} - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1}\cdot \mathbf{V}(\mathbf{n}, \mathbf{x}) ~. \end{align} $$ Also, taking the volume average of (13)$$_1$$, we have
 * $$ \text{(14)} \qquad

\langle \mathbf{P} \rangle = - \langle \left[\boldsymbol{S}(\mathbf{x}) - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} \rangle\cdot \mathbf{V}(\mathbf{n}, \mathbf{x}) ~. $$ Therefore, comparing (13)$$_2$$ and (14) and invoking the arbitrary nature $$\epsilon_0$$, we have

{  \left[\boldsymbol{S}_\text{eff} - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} = \langle \left[\boldsymbol{S}(\mathbf{x}) - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} \rangle ~. } $$ This relation provides us with a means of computing the effective permittivity of a layered medium oriented at a random angle (given by the normal $$\mathbf{n}$$) with respect to the coordinate basis.

Case 1: Simple or Rank-1 Laminate
Consider the Rank-1 laminate shown in Figure 3. The layers have permittivities alternating between $$\boldsymbol{\epsilon}_1$$ and $$\epsilon_2~\boldsymbol{\mathit{1}}$$. The volume fraction of phase $$1$$ is $$f_1$$ while that of phase $$2$$ is $$f_2$$ such that $$f_1 + f_2 = 1$$.

Recall that

\boldsymbol{S}(\mathbf{x}) = \epsilon_0~[\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}(\mathbf{x})]^{-1} ~. $$ Let us take the limit as $$\epsilon_0 \rightarrow \epsilon_2$$. Since $$\epsilon(\mathbf{x}) = \epsilon_2~\boldsymbol{\mathit{1}}$$ in phase $$2$$, we have

\boldsymbol{S}(\mathbf{x}) \rightarrow \infty \qquad \text{in phase}~2 ~. $$ Therefore,

\left[\boldsymbol{S}(\mathbf{x}) - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} \rightarrow 0 \qquad \text{in phase}~2 ~. $$ Hence, right hand side of

\left[\boldsymbol{S}_\text{eff} - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} = \langle \left[\boldsymbol{S}(\mathbf{x}) - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} \rangle $$ reduces to an average only over phase $$1$$. If we define

\boldsymbol{S}_1 := \epsilon_0~[\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_1]^{-1} \rightarrow \epsilon_2~[\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_1]^{-1} $$ we get
 * $$ \text{(15)} \qquad

\left[\boldsymbol{S}_\text{eff} - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} = f_1~\left[\boldsymbol{S}_1 - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right]^{-1} ~. $$ Taking the inverse of both sides of (15) gives

f_1~\left[\boldsymbol{S}_\text{eff} - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\right] = \boldsymbol{S}_1 - \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) $$ or,

f_1~\boldsymbol{S}_\text{eff} = \boldsymbol{S}_1 - (1 - f_1)~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) = \boldsymbol{S}_1 - f_2~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) ~. $$ Since

\boldsymbol{S}_\text{eff} = \epsilon_0~[\epsilon_0~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]^{-1} \rightarrow \epsilon_2~[\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]^{-1} $$ we then have

{  f_1~\epsilon_2~[\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]^{-1} = \epsilon_2~[\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_1]^{-1} - f_2~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) ~. } $$ This is the formula of Tartar, Murat, Lurie, and Cherkaev and can be shown to be equivalent to the Backus formula.


text:
Waves in composites and metamaterials/Hierarchical laminates and 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.

Hierarchical Laminates
In the previous lecture we found that, for rank-1 laminates, the effective permittivity can be calculated using the formula of Tartar-Murat-Lurie-Cherkaev. In this lecture we extend the ideas used to arrive at that formula to hierarchical laminates.

An example of a hierarchical laminate is shown in Figure 1. The idea of such materials goes back to Maxwell. In the rank-2 laminate shown in the figure there are two length scales which are assumed to be sufficiently separated so that the ideas in the previous lecture can be exploited. There has to be a separation of length scales so that the layer material can be replaced by its effective tensor.

Recall the Tartar-Murat-Lurie-Cherkaev formula for the effective permittivity of a rank-1 laminate:

f_1~[\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]^{-1} = [\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_1]^{-1} - \cfrac{f_2}{\epsilon_2}~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) ~. $$ By iterating this formula one gets, for a rank-$$m$$ laminate,
 * $$ \text{(1)} \qquad

f_1~[\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_\text{eff}]^{-1} = [\epsilon_2~\boldsymbol{\mathit{1}} - \boldsymbol{\epsilon}_1]^{-1} - \cfrac{f_2}{\epsilon_2}~\boldsymbol{M} $$ where

\boldsymbol{M} := \sum_{j=1}^m c_j~\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}_j) ~; c_j = \cfrac{f^{(j-1)} - f^{(j)}}{1 - f_1} \ge 0 $$ and $$m$$ is the number of laminates in the hierarchy, $$f^{(j)}$$ is the proportion of phase $$1$$ in a rank-$$j$$ laminate, and $$\mathbf{n}_j$$ is the orientation of the $$j$$-th laminate.

In particular,

f^{(m)} = f_1 ~;~ f^{(0)} = 1 \qquad \implies \qquad f^{(j-1)} > f^{(j)} ~. $$ Then,

\sum_{j=1}^m c_j = \cfrac{f^{(0)} - f^{(m)}}{1 - f_1} = 1~. $$

For a rank-3 laminate, if the normals $$\mathbf{n}_1$$, $$\mathbf{n}_2$$, and $$\mathbf{n}_3$$ are three orthogonal vectors, then

\boldsymbol{M} = c_1~\mathbf{n}_1\otimes\mathbf{n}_1 + c_2~\mathbf{n}_2\otimes\mathbf{n}_2 + c_3~\mathbf{n}_3\otimes\mathbf{n}_3 ~. $$ If we choose the $$f^{(j)}$$ s so that

c_1 = c_2 = c_3 = \frac{1}{3} $$ then

\boldsymbol{M} = \frac{1}{3}~(\mathbf{n}_1\otimes\mathbf{n}_1 + \mathbf{n}_2\otimes\mathbf{n}_2 +         \mathbf{n}_3\otimes\mathbf{n}_3) = \boldsymbol{\mathit{1}} ~. $$ In this case, equation (1) coincides with the solution for the Hashin sphere assemblage!

This implies that different geometries can have the same  $$\boldsymbol{\epsilon}_\text{eff}(\epsilon_2,\boldsymbol{\epsilon}_1)$$.

Is there a formula as simple as the Tartar-Murat-Lurie-Cherkaev formula when ε1 and ε2 are both anisotropic?
The answer is yes.

In this case we use an anisotropic reference material $$\boldsymbol{\epsilon}_0$$ and define the polarization as
 * $$ \text{(2)} \qquad

\mathbf{P}(\mathbf{x}) := [\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]\cdot\mathbf{E}(\mathbf{x}) = \mathbf{D}(\mathbf{x}) - \boldsymbol{\epsilon}_0\cdot\mathbf{E}(\mathbf{x}) ~. $$ The volume average of this field is given by

\langle \mathbf{P} \rangle := \langle \boldsymbol{\epsilon}\cdot\mathbf{E} \rangle - \boldsymbol{\epsilon}_0\cdot\langle \mathbf{E} \rangle = \langle \mathbf{D} \rangle - \boldsymbol{\epsilon}_0\cdot\langle \mathbf{E} \rangle ~. $$ Therefore, the difference between the field and its volume average is
 * $$ \text{(3)} \qquad

\mathbf{P}(\mathbf{x}) - \langle \mathbf{P} \rangle = [\mathbf{D}(\mathbf{x}) - \langle \mathbf{D} \rangle] - \boldsymbol{\epsilon}_0\cdot[\mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle] ~. $$ Let us introduce a new matrix $$\boldsymbol{\mathit{\Gamma}}(\mathbf{n})$$ defined through its action on a vector $$\mathbf{a}$$, i.e.,
 * $$ \text{(4)} \qquad

\mathbf{b} = \boldsymbol{\mathit{\Gamma}}(\mathbf{n}) \cdot \mathbf{a} \qquad\text{if and only if}\qquad \begin{cases} \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{b} = \mathbf{b} \\ \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot(\mathbf{a} - \boldsymbol{\epsilon}_0\cdot\mathbf{b}) = \boldsymbol{0} \\ \end{cases} $$ where $$\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) = \mathbf{n}\otimes\mathbf{n}$$ and projects parallel to $$\mathbf{n}$$. Therefore,

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{b} = \mathbf{b} \qquad \implies \qquad \alpha~\mathbf{n} = \mathbf{b} $$ where $$\alpha = \mathbf{b}\cdot\mathbf{n}$$. Also,

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot(\mathbf{a} - \boldsymbol{\epsilon}_0\cdot\mathbf{b}) = \boldsymbol{0} \qquad \implies \qquad \mathbf{a}\cdot\mathbf{n} - \alpha~(\mathbf{n}\cdot\boldsymbol{\epsilon}_0\cdot\mathbf{n}) = 0~. $$ Therefore,

\alpha = \cfrac{\mathbf{a}\cdot\mathbf{n}}{\mathbf{n}\cdot\boldsymbol{\epsilon}_0\cdot\mathbf{n}} \qquad \text{and} \qquad \mathbf{b} = \left(\cfrac{\mathbf{a}\cdot\mathbf{n}}{\mathbf{n}\cdot\boldsymbol{\epsilon}_0\cdot\mathbf{n}} \right)~\mathbf{n} = \cfrac{(\mathbf{n}\otimes\mathbf{n})\cdot\mathbf{a}}{\mathbf{n}\cdot\boldsymbol{\epsilon}_0\cdot\mathbf{n}} = \cfrac{\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{a}}{\mathbf{n}\cdot\boldsymbol{\epsilon}_0\cdot\mathbf{n}} ~. $$ From the definition of $$\boldsymbol{\mathit{\Gamma}}(\mathbf{n})$$ we then have
 * $$ \text{(5)} \qquad

\boldsymbol{\mathit{\Gamma}}(\mathbf{n}) = \cfrac{\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})}{\mathbf{n}\cdot\boldsymbol{\epsilon}_0\cdot\mathbf{n}} ~. $$ Taking the projection of both sides of equation (3) we get

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot[\mathbf{P}(\mathbf{x}) - \langle \mathbf{P} \rangle] = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot[\mathbf{D}(\mathbf{x}) - \langle \mathbf{D} \rangle] -    \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\{\boldsymbol{\epsilon}_0\cdot[\mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle]\} ~. $$ Now continuity of the normal component of $$\mathbf{D}$$ and the piecewise constant nature of the field implies that the normal component of $$\mathbf{D}$$ is constant. Therefore,

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot[\mathbf{D}(\mathbf{x}) - \langle \mathbf{D} \rangle] = [\mathbf{n}\cdot\mathbf{D}(\mathbf{x}) - \mathbf{n}\cdot\langle \mathbf{D} \rangle]~\mathbf{n} = \boldsymbol{0} ~. $$ Hence we have,
 * $$ \text{(6)} \qquad

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\left\{\mathbf{P}(\mathbf{x}) - \langle \mathbf{P} \rangle + \boldsymbol{\epsilon}_0\cdot[\mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle]\right\} = \boldsymbol{0} ~. $$ Recall from the previous lecture that
 * $$ \text{(7)} \qquad

\boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\mathbf{E}(\mathbf{x}) = \mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle + \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot\langle \mathbf{E} \rangle \quad \implies \quad \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n})\cdot[\mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle] = \mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle ~. $$ Since the conditions in (4) are satisfied with

\mathbf{a} := \mathbf{P}(\mathbf{x}) - \langle \mathbf{P} \rangle \quad \text{and} \quad \mathbf{b} := -[\mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle] $$ from the definition of $$\boldsymbol{\mathit{\Gamma}}$$ we then have
 * $$ \text{(8)} \qquad

-[\mathbf{E}(\mathbf{x}) - \langle \mathbf{E} \rangle] = \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\cdot[\mathbf{P}(\mathbf{x}) - \langle \mathbf{P} \rangle] ~. $$ Now, from equation (2) we have

\mathbf{E}(\mathbf{x}) = [\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]^{-1}\cdot\mathbf{P}(\mathbf{x}) ~. $$ Plugging this into (8) gives

-[\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]^{-1}\cdot\mathbf{P}(\mathbf{x}) + \langle \mathbf{E} \rangle = \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\cdot[\mathbf{P}(\mathbf{x}) - \langle \mathbf{P} \rangle] $$ or,

\{[\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}\cdot\mathbf{P}(\mathbf{x}) = \langle \mathbf{E} \rangle + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\cdot\langle \mathbf{P} \rangle ~. $$ Define

\mathbf{V} := \langle \mathbf{E} \rangle + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\cdot\langle \mathbf{P} \rangle $$ and note that this quantity is constant throughout the laminate. Therefore we can write

\{[\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}\cdot\mathbf{P}(\mathbf{x}) = \mathbf{V} $$ or

\mathbf{P}(\mathbf{x}) = \{[\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}^{-1}\cdot\mathbf{V} ~. $$ If we now take a volume average, we get
 * $$ \text{(9)} \qquad

\langle \mathbf{P} \rangle = \langle \{[\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]^{-1}+ \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}^{-1} \rangle \cdot\mathbf{V} ~. $$ Also, from the definition of $$\mathbf{P}(\mathbf{x})$$ we have

\langle \mathbf{P} \rangle = [\boldsymbol{\epsilon}_\text{eff} - \boldsymbol{\epsilon}_0]\cdot\langle \mathbf{E} \rangle \qquad \implies \qquad \langle \mathbf{E} \rangle = [\boldsymbol{\epsilon}_\text{eff} - \boldsymbol{\epsilon}_0]^{-1}\cdot\langle \mathbf{P} \rangle ~. $$ Therefore,

\mathbf{V} = \{[\boldsymbol{\epsilon}_\text{eff} - \boldsymbol{\epsilon}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}\cdot\langle \mathbf{P} \rangle $$ or,
 * $$ \text{(10)} \qquad

\langle \mathbf{P} \rangle = \{[\boldsymbol{\epsilon}_\text{eff} - \boldsymbol{\epsilon}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}^{-1}\cdot\mathbf{V}~. $$ Comparing equation (9) and (10) and invoking the arbitrariness of $$\boldsymbol{\epsilon}_0$$, we get
 * $$ \text{(11)} \qquad

{ \{[\boldsymbol{\epsilon}_\text{eff} - \boldsymbol{\epsilon}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}^{-1} = \langle \{[\boldsymbol{\epsilon}(\mathbf{x}) - \boldsymbol{\epsilon}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}^{-1} \rangle ~. } $$ This relation has a simple form and can be used when the phases are anisotropic.

For a simple (rank-1) laminate where $$\boldsymbol{\epsilon}_0 = \boldsymbol{\epsilon}_2$$, equation (11) reduces to

{ f_1~(\boldsymbol{\epsilon}_\text{eff} - \boldsymbol{\epsilon}_2)^{-1} = (\boldsymbol{\epsilon}_1 - \boldsymbol{\epsilon}_2)^{-1} + f_2~\boldsymbol{\mathit{\Gamma}}(\mathbf{n}) } $$ where

{ \boldsymbol{\mathit{\Gamma}}(\mathbf{n}) = \cfrac{\mathbf{n}\otimes\mathbf{n}}{\mathbf{n}\cdot\boldsymbol{\epsilon}_2\cdot\mathbf{n}} ~. } $$

Linear Elastic laminates
For elasticity, exactly the same analysis can be applied. In this case we introduce a reference stiffness tensor $$\boldsymbol{\mathsf{C}}_0$$ and define the second order polarization tensor as

\boldsymbol{P}(\mathbf{x}) = [\boldsymbol{\mathsf{C}}(\mathbf{x}) - \boldsymbol{\mathsf{C}}_0] : \boldsymbol{\varepsilon}(\mathbf{x}) $$ where the strain $$\boldsymbol{\varepsilon}(\mathbf{x})$$ is given by

\boldsymbol{\varepsilon}(\mathbf{x}) = \frac{1}{2}~[\boldsymbol{\nabla}\mathbf{u} + (\boldsymbol{\nabla}\mathbf{u})^T] ~. $$ Following the same process as before, we can show that the effective elastic stiffness of a hierarchical laminate can be determined from the formula
 * $$ \text{(12)} \qquad

{ \{[\boldsymbol{\mathsf{C}}_\text{eff} - \boldsymbol{\mathsf{C}}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}^{-1} = \langle \{[\boldsymbol{\mathsf{C}}(\mathbf{x}) - \boldsymbol{\mathsf{C}}_0]^{-1} + \boldsymbol{\mathit{\Gamma}}(\mathbf{n})\}^{-1} \rangle } $$ where (the components are in a rectangular Cartesian basis)

[\boldsymbol{\mathit{\Gamma}}(\mathbf{n})]_{ijlm} = \cfrac{1}{4}~\left[n_i~\{\boldsymbol{C}^{-1}(\mathbf{n})\}_{jl}~n_m + n_i~\{\boldsymbol{C}^{-1}(\mathbf{n})\}_{jm}~n_l + n_j~\{\boldsymbol{C}^{-1}(\mathbf{n})\}_{il}~n_m + n_j~\{\boldsymbol{C}^{-1}(\mathbf{n})\}_{im}~n_l\right] $$ and

\boldsymbol{C}^{-1}(\mathbf{n}) = \mathbf{n} \cdot \boldsymbol{\mathsf{C}}_0 \cdot \mathbf{n} ~. $$ Note that $$\boldsymbol{C}^{-1}(\mathbf{n})$$ has the same form as the acoustic tensor.

If $$\boldsymbol{\mathsf{C}}_0$$ is isotropic, i.e.,

\boldsymbol{\mathsf{C}}_0:\boldsymbol{\varepsilon} = \lambda_0~\text{tr}(\boldsymbol{\varepsilon})~\boldsymbol{\mathit{1}} + 2~\mu_0~\boldsymbol{\varepsilon} $$ where $$\lambda_0$$ is the Lame modulus and $$\mu_0$$ is the shear modulus, $$\boldsymbol{\mathit{\Gamma}}(\mathbf{n})$$ simplifies to

[\boldsymbol{\mathit{\Gamma}}(\mathbf{n})]_{ijlm} = \left(\cfrac{1}{\lambda_0 + 2~\mu_0} - \cfrac{1}{\mu_0}\right)~ n_i~n_j~n_l~n_m + \cfrac{1}{4~\mu_0} (n_i~\delta_{jl}~n_m + n_i~\delta_{jm}~n_l +    n_j~\delta_{il}~n_m + n_j~\delta_{im}~n_l) ~. $$

Hilbert Space Formulation
The methods discussed above can be generalized if we think in terms of a Hilbert space formalism. Recall that our goal is to find a general formula for $$\boldsymbol{\epsilon}_\text{eff}$$ and $$\boldsymbol{\mathsf{C}}_\text{eff}$$.

Let us consider a periodic material with unit cell $$Q$$. We will call such materials $$Q$$-periodic.

Electromagentism
Consider the Hilbert space $$\mathcal{H}$$ of square-integrable, $$Q$$-periodic, complex vector fields with the inner product

(\mathbf{a}, \mathbf{b}) = \int_Q \overline{\mathbf{a}(\mathbf{x})}\cdot\mathbf{b}(\mathbf{x})~\text{d}\mathbf{x} $$ where $$\mathbf{a}$$ and $$\mathbf{b}$$ are vector fields and $$\overline{(\bullet)}$$ denotes the complex conjugate. We can use Parseval's theorem to express the inner product in Fourier space as

(\mathbf{a}, \mathbf{b}) = \int_{\hat{Q}} \widehat{\mathbf{a}}(\mathbf{k})\cdot\widehat{\mathbf{b}}(\mathbf{k})~\text{d}\mathbf{k} $$ where $$\mathbf{k}$$ is the phase vector.

The Hilbert space $$\mathcal{H}$$ can be decomposed into three orthogonal subspaces.

 * 1) The subspace $$\mathcal{U}$$ of uniform fields, i.e., $$\mathbf{a}(\mathbf{x})$$ is independent of $$\mathbf{x}$$, or in Fourier space, $$\widehat{\mathbf{a}}(\mathbf{k}) = \boldsymbol{0}$$ unless $$\mathbf{k} = \boldsymbol{0}$$.
 * 2) The subspace $$\mathcal{J}$$ of zero divergence, zero average value fields, i.e., $$\boldsymbol{\nabla} \cdot \mathbf{a}(\mathbf{x}) = 0$$ and $$\langle \mathbf{a}(\mathbf{x}) \rangle = \boldsymbol{0}$$, or in Fourier space, $$\widehat{\mathbf{a}}(\boldsymbol{0}) = \boldsymbol{0}$$ and $$\mathbf{k}\cdot\widehat{\mathbf{a}}(\mathbf{k}) = 0$$.
 * 3) The subspace $$\mathcal{E}$$ of zero curl, zero average value fields, i.e., $$\boldsymbol{\nabla} \times \mathbf{a}(\mathbf{x}) = \boldsymbol{0}$$ and $$\langle \mathbf{a}(\mathbf{x}) \rangle = \boldsymbol{0}$$, or in Fourier space, $$\widehat{\mathbf{a}}(\boldsymbol{0}) = \boldsymbol{0}$$ and $$\widehat{\mathbf{a}}(\mathbf{k}) = \alpha(\mathbf{k})~\mathbf{k}$$.

Thus we can write

\mathcal{H} = \mathcal{U} \oplus \mathcal{J} \oplus \mathcal{E} ~. $$ In Fourier space, we can clearly see that

(\mathbf{a}, \mathbf{b}) = \boldsymbol{0} $$ if we choose $$\mathbf{a}$$ from any one of $$\mathcal{U}, \mathcal{J}, \mathcal{E}$$ and $$\mathbf{b}$$ from a a different subspace. Therefore the three subspaces are orthogonal.

Elasticity
Similarly, for elasticity, $$\mathcal{H}$$ is the Hilbert space of square-integrable, $$Q$$-periodic, complex valued, symmetric matrix valued fields with inner product

(\boldsymbol{A}, \boldsymbol{B}) = \int_Q \overline{\boldsymbol{A}(\mathbf{x})} : \boldsymbol{B}(\mathbf{x})~\text{d}\mathbf{x} ~. $$ In Fourier space, we have

(\boldsymbol{A}, \boldsymbol{B}) = \int_{\hat{Q}} \widehat{\boldsymbol{A}}(\mathbf{k}) : \widehat{\boldsymbol{B}}(\mathbf{k})~\text{d}\mathbf{k} ~. $$ Again we decompose the space $$\mathcal{H}$$ into three orthogonal subspaces $$\mathcal{U}$$, $$\mathcal{J}$$, and $$\mathcal{E}$$ where

 * 1) $$\mathcal{U}$$ is the subspace of uniform fields, i.e., $$\boldsymbol{A}(\mathbf{x})$$ is independent of $$\mathbf{x}$$, or in Fourier space, $$\widehat{\boldsymbol{A}}(\mathbf{k}) = \boldsymbol{\mathit{0}}$$ unless $$\mathbf{k} = \boldsymbol{0}$$.
 * 2) $$\mathcal{J}$$ is the subspace of zero divergence, zero average value fields, i.e., $$\boldsymbol{\nabla} \cdot \boldsymbol{A}(\mathbf{x}) = \boldsymbol{0}$$ and $$\langle \boldsymbol{A}(\mathbf{x}) \rangle = \boldsymbol{\mathit{0}}$$, or in Fourier space, $$\widehat{\boldsymbol{A}}(\boldsymbol{0}) = \boldsymbol{\mathit{0}}$$ and $$\widehat{\boldsymbol{A}}(\mathbf{k})\cdot\mathbf{k} = \boldsymbol{0}$$.
 * 3) $$\mathcal{E}$$ is the subspace of zero average "strain" fields, i.e., $$\langle \boldsymbol{A}(\mathbf{x}) \rangle = \boldsymbol{\mathit{0}}$$, or in Fourier space, $$\widehat{\boldsymbol{A}}(\boldsymbol{0}) = \boldsymbol{\mathit{0}}$$ and $$\widehat{\boldsymbol{A}}(\mathbf{k}) =  \mathbf{k}\otimes\mathbf{u}(\mathbf{k}) + \mathbf{u}(\mathbf{k})\otimes\mathbf{k}$$.

Problem of determining the effective tensor in an abstract setting
Let us first consider the problem of determining the effective permittivity. The approach will be to split relevant fields into components that belong to orthogonal subpaces of $$\mathcal{H}$$.

Since $$\boldsymbol{\nabla}\cdot\mathbf{D} = 0$$, we can split $$\mathbf{D}(\mathbf{x})$$ into two parts

\mathbf{D}(\mathbf{x}) = \mathbf{D}_0 + \mathbf{D}_1(\mathbf{x}) $$ where $$\mathbf{D}_0 \in \mathcal{U}$$ and $$\mathbf{D}_1 \in \mathcal{J}$$.

Also, since $$\boldsymbol{\nabla}\times\mathbf{E} = \boldsymbol{0}$$, we can split $$\mathbf{E}(\mathbf{x})$$ into two parts

\mathbf{E}(\mathbf{x}) = \mathbf{E}_0 + \mathbf{E}_1(\mathbf{x}) $$ where $$\mathbf{E}_0 \in \mathcal{U}$$ and $$\mathbf{E}_1 \in \mathcal{E}$$.

The constitutive relation linking $$\mathbf{D}$$ and $$\mathbf{E}$$ is

\mathbf{D}(\mathbf{x}) = \boldsymbol{\epsilon}(\mathbf{x})\cdot\mathbf{E}(\mathbf{x}) \equiv \boldsymbol{\epsilon}(\mathbf{E}) $$ where $$\boldsymbol{\epsilon}(\bullet)$$ can be thought of as an operator which is local in real space and maps $$\mathcal{H}$$ to $$\mathcal{H}$$. Therefore, we can write

\mathbf{D}_0 + \mathbf{D}_1 = \boldsymbol{\epsilon}(\mathbf{E}_0 + \mathbf{E}_1) ~. $$ The effective permittivity $$\boldsymbol{\epsilon}_\text{eff}$$ is defined through the relation

\mathbf{D}_0 = \boldsymbol{\epsilon}_\text{eff}(\mathbf{E}_0) ~. $$ Let $$\boldsymbol{\mathit{\Gamma}}_1$$ denote the projection operator that effects the projection of any vector in $$\mathcal{H}$$ onto the subspace $$\mathcal{E}$$. This projection is local in Fourier space. We can show that, if

\mathbf{b}(\mathbf{x}) = \boldsymbol{\mathit{\Gamma}}_1\cdot\mathbf{a}(\mathbf{x}) $$ then

\widehat{\mathbf{b}}(\mathbf{k}) = \boldsymbol{\mathit{\widehat{\Gamma}}}_1(\mathbf{k})\cdot\widehat{\mathbf{a}}(\mathbf{k}) $$ where

\boldsymbol{\mathit{\widehat{\Gamma}}}_1(\mathbf{k}) = \begin{cases} \cfrac{\mathbf{k}\otimes\mathbf{k}}{|\mathbf{k}|^2} = \boldsymbol{\mathit{\Gamma}}_1(\mathbf{n}) & \qquad\text{with} \mathbf{n} = \cfrac{\mathbf{k}}{|\mathbf{k}|} \text{if} \mathbf{k} \ne 0 \\ \boldsymbol{\mathit{0}} & \qquad\text{if}\mathbf{k} = 0 ~. \end{cases} $$ More generally, if we choose some reference matrix $$\boldsymbol{\epsilon}_0$$, we can define an operator $$\boldsymbol{\mathit{\Gamma}}$$ which is local in Fourier space via the relation

\mathbf{b}(\mathbf{x}) = \boldsymbol{\mathit{\Gamma}}\cdot\mathbf{a}(\mathbf{x}) $$ if and only if

\mathbf{b} \in \mathcal{E} \qquad \text{and} \qquad \mathbf{a} - \boldsymbol{\epsilon}_0(\mathbf{b}) \in \mathcal{U} \oplus \mathcal{J} \implies \boldsymbol{\nabla} \cdot (\mathbf{a} - \boldsymbol{\epsilon}_0(\mathbf{b})) = 0 ~. $$ In Fourier space,

\widehat{\mathbf{b}}(\mathbf{k}) = \boldsymbol{\mathit{\widehat{\Gamma}}}(\mathbf{k})\cdot\widehat{\mathbf{a}}(\mathbf{k}) $$ where

\boldsymbol{\mathit{\widehat{\Gamma}}}(\mathbf{k}) = \begin{cases} \cfrac{\mathbf{k}\otimes\mathbf{k}}{\mathbf{k}\cdot\boldsymbol{\epsilon}_0(\mathbf{k})} = \boldsymbol{\mathit{\Gamma}}(\mathbf{n}) & \qquad\text{with} \mathbf{n} = \cfrac{\mathbf{k}}{|\mathbf{k}|} \text{if} \mathbf{k} \ne 0 \\ \boldsymbol{\mathit{0}} & \qquad\text{if}\mathbf{k} = 0 ~. \end{cases} $$ In the next lecture we will derive relations for the effective tensors using these ideas.