# 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}}_{\perp }^{T}\cdot {\boldsymbol {R}}_{\perp })\cdot \langle \mathbf {D} \rangle =\lambda ~{\boldsymbol {R}}_{\perp }\cdot {\boldsymbol {\epsilon }}_{\text{eff}}^{-1}\cdot {\boldsymbol {R}}_{\perp }^{T}\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}}_{\perp }^{T}={\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.
* Figure 1. Polycrystal with randomly oriented lattice vectors.
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}}$
* Figure 2. Composite of two isotropic phases.
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{aligned}\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{aligned}}$
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.
* Figure 3. A laminate with direction of lamination $x_{1}$.
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 \neq \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 _{11}^{\text{eff}}}}~{\begin{bmatrix}-1&{\boldsymbol {\epsilon }}_{1t}^{\text{eff}}\\{\boldsymbol {\epsilon }}_{t1}^{\text{eff}}&~~\epsilon _{11}~{\boldsymbol {\epsilon }}_{tt}^{\text{eff}}-{\boldsymbol {\epsilon }}_{t1}^{\text{eff}}\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.

## Footnotes

* ↑ The following discussion is based on Milton02.