# Waves in composites and metamaterials/Bloch waves and the quasistatic limit

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.

## Bloch Theorem

In the previous lecture we showed that Maxwell's equations at fixed frequency can be formulated in terms of the fields $\mathbf {D}$ and $\mathbf {B}$ as
    ${\text{(1)}}\qquad {\boldsymbol {\nabla }}\cdot \mathbf {D} =0~;~~{\boldsymbol {\nabla }}\cdot \mathbf {B} =0~;~~i~{\boldsymbol {\nabla }}\times \left({\cfrac {\mathbf {D} }{\epsilon }}\right)=\omega ~\mathbf {B} ~;~~-i~{\boldsymbol {\nabla }}\times \left({\cfrac {\mathbf {B} }{\mu }}\right)=\omega ~\mathbf {D} ~.$
Equations (1) suggest that we should look for solutions $\mathbf {D}$ and $\mathbf {B}$ in the space of divergence-free fields such that
    ${\text{(2)}}\qquad {\mathcal {L}}{\begin{bmatrix}\mathbf {D} \\\mathbf {B} \end{bmatrix}}=\omega ~{\begin{bmatrix}\mathbf {D} \\\mathbf {B} \end{bmatrix}}\qquad {\text{or}}\qquad {({\mathcal {L}}-\omega ~{\boldsymbol {1}}){\begin{bmatrix}\mathbf {D} \\\mathbf {B} \end{bmatrix}}={\boldsymbol {0}}}$
where the operator ${\mathcal {L}}$ is given by
    ${\text{(3)}}\qquad {{\mathcal {L}}:={\begin{bmatrix}0&-i~{\boldsymbol {\nabla }}\times \mu ^{-1}\\i~{\boldsymbol {\nabla }}\times \epsilon ^{-1}&0\end{bmatrix}}~.}$
If the medium is such that the permittivity $\epsilon (\mathbf {x} )$ and the permeability $\mu (\mathbf {x} )$ are periodic, i.e.,
    $\epsilon (\mathbf {x} )=\epsilon (\mathbf {x} +\mathbf {R} )~;~~\mu (\mathbf {x} )=\mu (\mathbf {x} +\mathbf {R} )$
where $\mathbf {R}$ is a lattice vector (see Figure 1) then the operator ${\mathcal {L}}$ has the same periodicity as the medium.
* Figure 1. Lattice vector in a periodic medium.
Also recall the translation operator ${\mathcal {T}}_{R}$ defined as
    ${{\mathcal {T}}_{R}{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}={\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} )\\\mathbf {B} (\mathbf {x} +\mathbf {R} )\end{bmatrix}}~.}$
Periodicity of the medium implies that ${\mathcal {T}}_{R}$ commutes with ${\mathcal {L}}$, i.e.,
    ${{\mathcal {T}}_{R}~{\mathcal {L}}={\mathcal {L}}~{\mathcal {T}}_{R}~.}$
The translation operator is unitary, i.e.,
    ${{\mathcal {T}}_{R}~{\mathcal {T}}_{R}^{T}={\mathcal {T}}_{R}~{\mathcal {T}}_{-R}={\boldsymbol {1}}~.}$
This means that the adjoint operator ${\mathcal {T}}_{R}^{T}$ is equal to the inverse operator ${\mathcal {T}}_{-R}$.
The translation operator also commutes, i.e.,
    ${{\mathcal {T}}_{R}~{\mathcal {T}}_{R'}={\mathcal {T}}_{R'}~{\mathcal {T}}_{R}={\mathcal {T}}_{R+R'}~.}$
Also, since ${\mathcal {L}}$ and ${\mathcal {T}}_{R}$ commute, the operators ${\mathcal {L}}-\omega {\boldsymbol {1}}$ and ${\mathcal {T}}_{R}$ must also commute. This implies that
    ${\mathcal {T}}_{R}~({\mathcal {L}}-\omega ~{\boldsymbol {1}})~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}={\boldsymbol {0}}=({\mathcal {L}}-\omega ~{\boldsymbol {1}})~{\mathcal {T}}_{R}~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}=({\mathcal {L}}-\omega ~{\boldsymbol {1}})~~{\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} )\\\mathbf {B} (\mathbf {x} +\mathbf {R} )\end{bmatrix}}~.$
Hence the eigenstates of ${\mathcal {L}}-\omega ~{\boldsymbol {1}}$ and the eigenstates of ${\mathcal {T}}_{R}$ lie in the same space. Therefore, any solution can be expressed in fields which are simultaneously eigenstates of all the ${\mathcal {T}}_{R}$, i.e., these eigenstates have the property
    ${\text{(4)}}\qquad {[{\mathcal {T}}_{R}-c(\mathbf {R} )~{\boldsymbol {1}}]{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}={\boldsymbol {0}}~.}$
Since ${\mathcal {T}}_{R}~{\mathcal {T}}_{R'}={\mathcal {T}}_{R+R'}$, we have
    $c(\mathbf {R} )~c(\mathbf {R} ')=c(\mathbf {R} +\mathbf {R} ')~.$
So it suffices to know $c(\mathbf {R} )$ when $\mathbf {R} =\mathbf {a} _{1},\mathbf {a} _{2},\mathbf {a} _{3}$ where the $\mathbf {a} _{i}$'s are the primitive vectors of the lattice, i.e.,
    $\epsilon (\mathbf {x} +\mathbf {a} _{j})=\epsilon (\mathbf {x} )~;~~\mu (\mathbf {x} +\mathbf {a} _{j})=\mu (\mathbf {x} )~.$
Let us assume that
    $c(\mathbf {a} _{j})=e^{2\pi ~i~\alpha _{j}}\qquad j=1,2,3$
for a suitable choice of $\alpha _{j}$.
Then for any lattice vector
    $\mathbf {R} =n_{1}~\mathbf {a} _{1}+n_{2}~\mathbf {a} _{2}+n_{3}~\mathbf {a} _{3}$
we have
    ${\begin{aligned}c(\mathbf {R} )&=c(n_{1}~\mathbf {a} _{1}+n_{2}~\mathbf {a} _{2}+n_{3}~\mathbf {a} _{3})=c(n_{1}~\mathbf {a} _{1})~c(n_{2}~\mathbf {a} _{2})~c(n_{3}~\mathbf {a} _{3})\\&=c\left(\sum _{i=1}^{n_{1}}\mathbf {a} _{1}\right)~c\left(\sum _{i=1}^{n_{2}}\mathbf {a} _{2}\right)~c\left(\sum _{i=1}^{n_{3}}\mathbf {a} _{3}\right)=[c(\mathbf {a} _{1})]^{n_{1}}~[c(\mathbf {a} _{2})]^{n_{2}}~[c(\mathbf {a} _{3})]^{n_{3}}\\&=e^{2\pi i\alpha _{1}~n_{1}}~e^{2\pi i\alpha _{2}~n_{2}}~e^{2\pi i\alpha _{3}~n_{3}}=e^{2\pi i(\alpha _{1}~n_{1}+\alpha _{2}~n_{2}+\alpha _{3}~n_{3})}\end{aligned}}$
or,
    $c(\mathbf {R} )=e^{2\pi i(\alpha _{1}~n_{1}+\alpha _{2}~n_{2}+\alpha _{3}~n_{3})}~.$
Define a vector
    $\mathbf {k} :=\alpha _{1}~\mathbf {b} _{1}+\alpha _{2}~\mathbf {b} _{2}+\alpha _{3}~\mathbf {b} _{3}$
where the vectors $\mathbf {b} _{i}$ are the reciprocal lattice vectors satisfying
    $\mathbf {b} _{i}\cdot \mathbf {a} _{j}=2\pi ~\delta _{ij}~.$
Then,
    $\mathbf {k} \cdot \mathbf {R} =(\alpha _{1}~\mathbf {b} _{1}+\alpha _{2}~\mathbf {b} _{2}+\alpha _{3}~\mathbf {b} _{3})\cdot (n_{1}~\mathbf {a} _{1}+n_{2}~\mathbf {a} _{2}+n_{3}~\mathbf {a} _{3})=\alpha _{1}~n_{1}~2~\pi +\alpha _{2}~n_{2}~2~\pi +\alpha _{3}~n_{3}~2~\pi$
or,
    $\mathbf {k} \cdot \mathbf {R} =2\pi (\alpha _{1}~n_{1}+\alpha _{2}~n_{2}+\alpha _{3}~n_{3})~.$
Therefore, we have
    $c(\mathbf {R} )=e^{i\mathbf {k} \cdot \mathbf {R} }~.$
Plugging this expression into (4), we get
    $[{\mathcal {T}}_{R}-e^{i\mathbf {k} \cdot \mathbf {R} }~{\boldsymbol {1}}]{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}={\boldsymbol {0}}$
or,
    ${\text{(5)}}\qquad {{\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} )\\\mathbf {B} (\mathbf {x} +\mathbf {R} )\end{bmatrix}}=e^{i\mathbf {k} \cdot \mathbf {R} }~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}~.}$
Equation (5) is called the Bloch condition.
In summary, the solutions to the electromagnetic equations in a periodic medium can be expressed in Bloch waves where each Bloch wave is a time harmonic solution to the electromagnetic equations which in addition satisfies the Bloch condition for all lattice vectors $\mathbf {R}$ and for some appropriate choice of $\mathbf {k}$.
Note that for any vector $\mathbf {x}$, the Bloch condition implies that
    $e^{-i\mathbf {k} \cdot (\mathbf {x} +\mathbf {R} )}{\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} )\\\mathbf {B} (\mathbf {x} +\mathbf {R} )\end{bmatrix}}=e^{i\mathbf {k} \cdot \mathbf {R} -i\mathbf {k} \cdot \mathbf {R} -i\mathbf {k} \cdot \mathbf {x} }~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}=e^{-i\mathbf {k} \cdot \mathbf {x} }~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}~.$
Therefore the quantity
    $e^{-i\mathbf {k} \cdot \mathbf {x} }~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}$
is periodic.

## Quasistatic Limit

Let us now consider the solution of Maxwell's equation in periodic media in the quasistatic limit. Consider the periodic medium shown in Figure 2. The lattice spacing is $\eta$.
* Figure 2. Periodic medium with $x$ and $y$ spaces.
Define
    $\epsilon _{\eta }(\mathbf {x} ):=\epsilon \left({\cfrac {\mathbf {x} }{\eta }}\right)=\epsilon (\mathbf {y} )~;~~\mu _{\eta }(\mathbf {x} ):=\mu \left({\cfrac {\mathbf {x} }{\eta }}\right)=\mu (\mathbf {y} )~.$
These are periodic functions, i.e.,
    $\epsilon (\mathbf {y} +\mathbf {a} _{i})=\epsilon (\mathbf {y} )~;~~\mu (\mathbf {y} +\mathbf {a} _{i})=\mu (\mathbf {y} )$
where $\mathbf {a} _{i}$ are the primitive lattice vectors. We may also write these periodicity conditions as
    $\epsilon _{\eta }(\mathbf {x} +\eta \mathbf {a} _{i})=\epsilon _{\eta }(\mathbf {x} )~;~~\mu _{\eta }(\mathbf {x} +\eta \mathbf {a} _{i})=\mu _{\eta }(\mathbf {x} )~.$
Similarly, define
    $\mathbf {D} _{\eta }(\mathbf {x} ):=\mathbf {D} (\mathbf {y} )~;~~\mathbf {B} _{\eta }(\mathbf {x} ):=\mathbf {B} (\mathbf {y} )~;~~\mathbf {E} _{\eta }(\mathbf {x} ):=\mathbf {E} (\mathbf {y} )~;~~\mathbf {H} _{\eta }(\mathbf {x} ):=\mathbf {H} (\mathbf {y} )~.$
Then Maxwell's equations can be written as
    ${\text{(6)}}\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}}~.$
Let us look for Bloch wave solutions of the form
    $\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} )$
where $\mathbf {e} _{\eta },\mathbf {d} _{\eta },\mathbf {h} _{\eta },\mathbf {b} _{\eta }$ have the same periodicity as $\epsilon$ and $\mu$, i.e.,
    $\mathbf {e} _{\eta }(\mathbf {x} +\eta ~\mathbf {a} _{i})=\mathbf {e} (\mathbf {x} )~;~~\mathbf {d} _{\eta }(\mathbf {x} +\eta ~\mathbf {a} _{i})=\mathbf {d} (\mathbf {x} )~;~~\mathbf {h} _{\eta }(\mathbf {x} +\eta ~\mathbf {a} _{i})=\mathbf {h} (\mathbf {x} )~;~~\mathbf {b} _{\eta }(\mathbf {x} +\eta ~\mathbf {a} _{i})=\mathbf {b} (\mathbf {x} )~.$
From the constitutive relations, we get
    $\mathbf {d} _{\eta }(\mathbf {x} )=\epsilon _{\eta }(\mathbf {x} )~\mathbf {e} _{\eta }(\mathbf {x} )~;~~\mathbf {b} _{\eta }(\mathbf {x} )=\mu _{\eta }(\mathbf {x} )~\mathbf {h} _{\eta }(\mathbf {x} )~.$
Recall that, for periodic media, Maxwell's equations may be expressed as
    $({\mathcal {L}}-\omega ~{\boldsymbol {1}}){\begin{bmatrix}\mathbf {D} \\\mathbf {B} \end{bmatrix}}={\boldsymbol {0}}~.$
Here $\omega$ is an eigenvalue of ${\mathcal {L}}$. However, ${\mathcal {L}}$ depends on $\omega$ via $\epsilon (\omega )$ and $\mu (\omega )$. {\bf Bloch wave solutions do not exists unless $\omega$ takes one of a discrete set of values.}
Let these discrete values be
    $\omega =\omega _{\eta }^{j}(\mathbf {k} )$
where the superscript $j$ labels the solution branches.
Let us see what the Bloch wave solutions reduce to as $\eta \rightarrow 0$. Following standard multiple scale analysis, let us assume that the periodic complex fields have the expansions
    ${\text{(7)}}\qquad {\begin{aligned}\mathbf {e} _{\eta }(\mathbf {x} )&=\mathbf {e} _{0}(\mathbf {y} )+\eta ~\mathbf {a} _{1}(\mathbf {y} )+\eta ^{2}~\mathbf {a} _{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{aligned}}$
Let us also assume that the dependence of $\omega$ on $\eta$ and $\mathbf {k}$ has an expansion of the form
    ${\text{(8)}}\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$
Plugging (8) and (7) into (6) gives
    ${\text{(9)}}\qquad {\begin{aligned}{\boldsymbol {\nabla }}\cdot \left(e^{i\eta ~\mathbf {k} \cdot \mathbf {y} }~[\mathbf {d} _{0}(\mathbf {y} )+\eta ~\mathbf {d} _{1}(\mathbf {y} )+\dots ]\right)&=0\\{\boldsymbol {\nabla }}\cdot \left(e^{i\eta ~\mathbf {k} \cdot \mathbf {y} }~[\mathbf {b} _{0}(\mathbf {y} )+\eta ~\mathbf {b} _{1}(\mathbf {y} )+\dots ]\right)&=0\\{\boldsymbol {\nabla }}\times \left(e^{i\eta ~\mathbf {k} \cdot \mathbf {y} }~[\mathbf {e} _{0}(\mathbf {y} )+\eta ~\mathbf {e} _{1}(\mathbf {y} )+\dots ]\right)-i[\omega _{0}^{j}+\eta ~\omega _{1}^{j}+\dots ][\mathbf {b} _{0}(\mathbf {y} )+\eta ~\mathbf {b} _{1}(\mathbf {y} )+\dots ]&={\boldsymbol {0}}\\{\boldsymbol {\nabla }}\times \left(e^{i\eta ~\mathbf {k} \cdot \mathbf {y} }~[\mathbf {h} _{0}(\mathbf {y} )+\eta ~\mathbf {h} _{1}(\mathbf {y} )+\dots ]\right)+i[\omega _{0}^{j}+\eta ~\omega _{1}^{j}+\dots ][\mathbf {d} _{0}(\mathbf {y} )+\eta ~\mathbf {d} _{1}(\mathbf {y} )+\dots ]&={\boldsymbol {0}}~.\end{aligned}}$
Define
    ${\boldsymbol {\nabla }}_{y}:=\left({\frac {\partial }{\partial y_{1}}},{\frac {\partial }{\partial y_{2}}},{\frac {\partial }{\partial y_{3}}}\right)~.$
Then, for a vector field $\mathbf {v} (\mathbf {y} )$, using the chain rule we get
    ${\text{(10)}}\qquad {\boldsymbol {\nabla }}\cdot \mathbf {v} (\mathbf {y} )={\cfrac {1}{\eta }}~{\boldsymbol {\nabla }}_{y}\cdot \mathbf {v} (\mathbf {y} )~;~~{\boldsymbol {\nabla }}\times \mathbf {v} (\mathbf {y} )={\cfrac {1}{\eta }}~{\boldsymbol {\nabla }}_{y}\times {\mathbf {v} (\mathbf {y} )}~.$
Using definitions (10) in (9) and collecting terms of order $1/\eta$ gives
    ${\text{(11)}}\qquad {\begin{aligned}{\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}}\end{aligned}}$
These are the solutions in the quasistatic limit. Also, from the constitutive equations
    $\mathbf {d} _{0}(\mathbf {y} )=\epsilon (\mathbf {y} )~\mathbf {e} _{0}(\mathbf {y} )~;~~\mathbf {b} _{0}(\mathbf {y} )=\mu (\mathbf {y} )~\mathbf {h} _{0}(\mathbf {y} )~.$
Similarly, collecting terms of order 1 from the expanded Maxwell's equations (9) we get
    ${\text{(12)}}\qquad {\begin{aligned}i~\mathbf {k} \cdot \mathbf {d} _{0}(\mathbf {y} )+{\boldsymbol {\nabla }}_{y}\cdot \mathbf {d} _{1}(\mathbf {y} )=0\\i~\mathbf {k} \cdot \mathbf {b} _{0}(\mathbf {y} )+{\boldsymbol {\nabla }}_{y}\cdot \mathbf {b} _{1}(\mathbf {y} )=0\\i~\mathbf {k} \times \mathbf {e} _{0}(\mathbf {y} )+{\boldsymbol {\nabla }}_{y}\times {\mathbf {e} _{1}(\mathbf {y} )}-i~\omega _{0}^{j}~\mathbf {b} _{0}(\mathbf {y} )={\boldsymbol {0}}\\i~\mathbf {k} \times \mathbf {h} _{0}(\mathbf {y} )+{\boldsymbol {\nabla }}_{y}\times {\mathbf {h} _{1}(\mathbf {y} )}+i~\omega _{0}^{j}~\mathbf {d} _{0}(\mathbf {y} )={\boldsymbol {0}}~.\end{aligned}}$
Since $\mathbf {d} _{1}(\mathbf {y} ),\mathbf {b} _{1}(\mathbf {y} ),\mathbf {e} _{1}(\mathbf {y} ),\mathbf {h} _{1}(\mathbf {y} )$ are periodic, this implies that
    ${\text{(13)}}\qquad {\begin{aligned}\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}}~.\end{aligned}}$
where $\langle \bullet \rangle$ is the volume average over the unit cell. So a necessary condition that equations (12) have a solution is that
    ${\text{(14)}}\qquad {\begin{aligned}i~\mathbf {k} \cdot \langle \mathbf {d} _{0}(\mathbf {y} )\rangle =0\\i~\mathbf {k} \cdot \langle \mathbf {b} _{0}(\mathbf {y} )\rangle =0\\i~\mathbf {k} \times \langle \mathbf {e} _{0}(\mathbf {y} )\rangle -i~\omega _{0}^{j}~\langle \mathbf {b} _{0}(\mathbf {y} )\rangle ={\boldsymbol {0}}\\i~\mathbf {k} \times \langle \mathbf {h} _{0}(\mathbf {y} )\rangle +i~\omega _{0}^{j}~\langle \mathbf {d} _{0}(\mathbf {y} )\rangle ={\boldsymbol {0}}~.\end{aligned}}$
Note that the second pair of (14) implies the first pair.

## Footnotes

* ↑ The following discussion is based on Ashcroft76 (p. 133-139).
* ↑ We can see that the two operators commute by working out the operations. Thus, ${\mathcal {T}}_{R}~{\mathcal {L}}~{\begin{bmatrix}\mathbf {D} \\\mathbf {B} \end{bmatrix}}={\mathcal {T}}_{R}~{\begin{bmatrix}-i~{\boldsymbol {\nabla }}\times [\mu ^{-1}(\mathbf {x} )~\mathbf {B} (\mathbf {x} )]\\i~{\boldsymbol {\nabla }}\times [\epsilon ^{-1}(\mathbf {x} )~\mathbf {D} (\mathbf {x} )]\end{bmatrix}}={\begin{bmatrix}-i~{\boldsymbol {\nabla }}\times [\mu ^{-1}(\mathbf {x} +\mathbf {R} )~\mathbf {B} (\mathbf {x} +\mathbf {R} )]\\i~{\boldsymbol {\nabla }}\times [\epsilon ^{-1}(\mathbf {x} +\mathbf {R} )~\mathbf {D} (\mathbf {x} +\mathbf {R} )]\end{bmatrix}}={\mathcal {L}}~{\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} )\\\mathbf {B} (\mathbf {x} +\mathbf {R} )\end{bmatrix}}={\mathcal {L}}~{\mathcal {T}}_{R}~{\begin{bmatrix}\mathbf {D} \\\mathbf {B} \end{bmatrix}}~.$
* ↑ We can see that the translation operator commutes by working out the operations. Thus, ${\mathcal {T}}_{R}~{\mathcal {T}}_{R'}~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}={\mathcal {T}}_{R}~{\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} ')\\\mathbf {B} (\mathbf {x} +\mathbf {R} ')\end{bmatrix}}={\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} '+\mathbf {R} )\\\mathbf {B} (\mathbf {x} +\mathbf {R} '+\mathbf {R} )\end{bmatrix}}={\mathcal {T}}_{R+R'}~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}={\mathcal {T}}_{R'}~{\begin{bmatrix}\mathbf {D} (\mathbf {x} +\mathbf {R} )\\\mathbf {B} (\mathbf {x} +\mathbf {R} )\end{bmatrix}}={\mathcal {T}}_{R'}~{\mathcal {T}}_{R}~{\begin{bmatrix}\mathbf {D} (\mathbf {x} )\\\mathbf {B} (\mathbf {x} )\end{bmatrix}}~.$
* ↑ The following discussion is based on Milton02