Waves in composites and metamaterials/Elastodynamics and electrodynamics

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.

Dissipation

Recall from the previous lecture that the average rate of work done in a cycle of oscillation of material with frequency dependent mass is

$\begin{array}{cc}𝒫& =\frac{\omega }{2\pi }\underset{0}{\overset{2\pi /\omega }{\int }}𝐅(t)\cdot 𝐕(t)\text{d}t\\ & =\frac{\text{Re}(\widehat{𝐅})\cdot \text{Re}(\widehat{𝐕})+\text{Im}(\widehat{𝐅})\cdot \text{Im}(\widehat{𝐕})}{2}\\ & =\omega [\text{Re}(\widehat{𝐕})\cdot \text{Im}[𝑴(\omega )]\cdot \text{Re}(\widehat{𝐕})+\text{Im}(\widehat{𝐕})\cdot \text{Im}[𝑴(\omega )]\cdot \text{Im}(\widehat{𝐕})].\end{array}$

This quadratic form will be non-negative for all choices of $\widehat{𝐕}$ if and only if $\text{Im}(𝑴(\omega ))$ is positive semidefinite for all real $\omega >0$. Therefore, a restriction on the behavior of such materials is that

$\text{Im}[𝑴(\omega )]\ge 0.$

Similarly, for electrodynamics, the average power dissipated into heat is given by

$\text{(1)}𝒫=\frac{\omega }{2\pi }\underset{0}{\overset{2\pi /\omega }{\int }}[𝐄(t)\cdot \frac{\partial 𝐃(t)}{\partial t}+𝐇(t)\cdot \frac{\partial 𝐁(t)}{\partial t}]dt.$

In this case, the quantity $𝐄$ is equivalent to the voltage and the quantity rate of change of electrical displacement $\partial 𝐃/\partial t$ is equivalent to the current (recall that in electrostatics the power is given by $𝒫=VI$). In addition, we also have a contribution due to magnetic induction.

Let us assume that the fields can be expressed in harmonic form, i.e.,

$𝐄(t)=\text{Re}[\widehat{𝐄}{e}^{-i\omega t}];𝐇(t)=\text{Re}[\widehat{𝐇}{e}^{-i\omega t}]$

or equivalently as

$𝐄(t)=\text{Re}(\widehat{𝐄})\cos (\omega t)+\text{Im}(\widehat{𝐄})\sin (\omega t);𝐇(t)=\text{Re}(\widehat{𝐇})\cos (\omega t)+\text{Im}(\widehat{𝐇})\sin (\omega t).$

Also, recall that,

$\frac{\partial 𝐃}{\partial t}=\mathrm{\nabla }\times 𝐇=-i\omega \mathit{ϵ}\cdot 𝐄(t)\text{and}\frac{\partial 𝐁}{\partial t}=-\mathrm{\nabla }\times 𝐄=-i\omega \mathit{\mu }\cdot 𝐇(t).$

Therefore, for real $\omega$ and real $𝒫$, we can write equation (1) as (with the substitution $z=\omega t$),

$\begin{array}{cc}𝒫=\frac{\omega }{2\pi }\underset{0}{\overset{2\pi }{\int }}& [\text{Re}(\widehat{𝐄})\cos z+\text{Im}(\widehat{𝐄})\sin z]\cdot \{\text{Im}(\mathit{ϵ})\cdot [\text{Re}(\widehat{𝐄})\cos z+\text{Im}(\widehat{𝐄})\sin z]\}+\\ & [\text{Re}(\widehat{𝐇})\cos z+\text{Im}(\widehat{𝐇})\sin z]\cdot \{\text{Im}(\mathit{\mu })\cdot [\text{Re}(\widehat{𝐇})\cos z+\text{Im}(\widehat{𝐇})\sin z]\}\text{d}z\end{array}$

Expanding out, and using the fact that

$\underset{0}{\overset{2\pi }{\int }}{\cos }^{2}z\text{d}z=\underset{0}{\overset{2\pi }{\int }}{\sin }^{2}z\text{d}z=\pi \text{and}\underset{0}{\overset{2\pi }{\int }}\sin (2z)\text{d}z=0$

we have,

$\text{(2)}𝒫=\frac{\omega }{2}[\text{Re}(\widehat{𝐄})\cdot \text{Im}(\mathit{ϵ})\cdot \text{Re}(\widehat{𝐄})+\text{Im}(\widehat{𝐄})\cdot \text{Im}(\mathit{ϵ})\cdot \text{Im}(\widehat{𝐄})+\text{Re}(\widehat{𝐇})\cdot \text{Im}(\mathit{\mu })\cdot \text{Re}(\widehat{𝐇})+\text{Im}(\widehat{𝐇})\cdot \text{Im}(\mathit{\mu })\cdot \text{Im}(\widehat{𝐇})].$

Since $\omega \ge 0$ and the power $𝒫\ge 0$, the quadratic forms in equation (2) require that

$\text{Im}(\mathit{ϵ})>0\text{and}\text{Im}(\mathit{\mu })>0.$

Note that if the permittivity is expressed as

$\mathit{ϵ}={\mathit{ϵ}}_0+i\frac{\mathit{\sigma }}{\omega }$

the requirement $\text{Im}(\mathit{ϵ})>0$ implies that the conductivity $\mathit{\sigma }>0$. Therefore, if the conductivity is greater than zero, there will be dissipation.

Brief introduction to elastodynamics

A concise introduction to the theory of elasticity can be found in Atkin80. In this section, we consider the linear theory of elasticity for infinitesimal strains and small displacements.

Consider the body ($\mathrm{\Omega }$) shown in Figure~1. Let $\mathrm{\Gamma }$ be a subpart of the body (in the interior of $\mathrm{\Omega }$ or sharing a part of the surface of $\mathrm{\Omega }$). Postulate the existence of a force $𝐭(𝐱,𝐧)$ per unit area on the surface of $\mathrm{\Gamma }$ where $𝐧$ is the outward unit normal to the surface of $\mathrm{\Gamma }$. Then $𝐭$ is the force exerted on $\mathrm{\Gamma }$ by the material outside $\mathrm{\Gamma }$ or by surface tractions.

File:PotatoStress.jpg Figure 1. Illustration of the concept of stress.

From the balance of forces on a small tetrahedron ($\mathrm{\Gamma }$), we can show that $𝐭(𝐱,𝐧)$ is linear in $𝐧$. Therefore,

$𝐭=\mathit{\sigma }\cdot 𝐧$

where $\mathit{\sigma }$ is a second-order tensor called the stress tensor.

Since the tetrahedron cannot rotate at infinite velocity as its size goes to zero (conservation of angular momentum), we can show that the stress tensor is symmetric, i.e.,

$\mathit{\sigma }={\mathit{\sigma }}^T.$

In particular, for a fluid,

$\mathit{\sigma }=-p1$

where $p$ is the pressure.

Let us assume that the stress depends only on the strain (and not on strain gradients or strain rates), where the strain is defined as

$\text{(3)}\mathit{ϵ}=\frac{1}{2}[\mathrm{\nabla }𝐮+(\mathrm{\nabla }𝐮{)}^T].$

Here $𝐮(𝐱,t)$ is the displacement field. Note that a gradient of the displacement field is used to define the strain because rigid body motions should not affect $\mathit{\sigma }$ and a rigid body rotation gives zero strains (for small displacements).

Assume that $\mathit{\sigma }$ depends linearly on $\mathit{ϵ}$ so that

$\mathit{\sigma }(𝐱,t)=\int \text{d}{𝐱}^{\prime }[\int {𝖪}_{\epsilon }(𝐱,{𝐱}^{\prime },t^{\prime }-t):\mathit{ϵ}(𝐱,t^{\prime })\text{d}{t}^{\prime }].$

Note that this assumption ignores preexisting internal stresses such as those found in prestressed concrete. If the material can be approximated as being local, then

$\text{(4)}\mathit{\sigma }(𝐱,t)=\int {𝖪}_{\epsilon }(𝐱,t^{\prime }-t):\mathit{ϵ}(𝐱,t^{\prime })\text{d}{t}^{\prime }.$

Taking the Fourier transforms of equation (4), we get

$\text{(5)}\hat{\mathit{\sigma }}(𝐱,\omega )=𝖢(𝐱,\omega ):\widehat{\mathit{ϵ}}(𝐱,\omega )$

where

$𝖢(𝐱,\omega )=\int {𝖪}_{\epsilon }(𝐱,\tau )e^{-i\omega \tau }\text{d}\tau ;\tau =t^{\prime }-t.$

In index notation, equation (5) can be written as

${\hat{\sigma }}_{ij}=C_{ijkl}{\widehat{ϵ}}_{kl}.$

Causality implies that stresses at time $t$ can only depend on strains of previous times, i.e., if $t^{\prime }\le t$ or $\tau \le 0$. Therefore,

${𝖪}_{\epsilon }(𝐱,\tau )=\mathrm{𝟢}\text{if}\tau >0.$

This in turn implies that the integral converges only if $\text{Im}(\omega )>0$, i.e., $𝖢(𝐱,\omega )$ is analytic when $\text{Im}(\omega )>0$.

In the absence of body forces, the equation of motion of the body can be written as

$\text{(6)}\mathrm{\nabla }\cdot \mathit{\sigma }=\rho (𝐱)\frac{{\partial }^2𝐮}{\partial t^2}$

where $\rho$ is the mass density, $\mathrm{\nabla }\cdot \mathit{\sigma }$ is the internal force per unit volume, and ${\partial }^2𝐮/\partial t^2$ is the acceleration. Hence, this is just the expression of Newton's second law for continuous systems.

For a material which has a frequency dependent mass, equation (6) may be written as

$\text{(7)}\mathrm{\nabla }\cdot \mathit{\sigma }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\rho (𝐱,t-t^{\prime })\frac{{\partial }^2𝐮}{\partial t^2}\text{d}{t}^{\prime }$

where causality implies that if $t^{\prime }>t$ then $\rho =0$.

Taking the Fourier transform of equation (7), we get

$\text{(8)}\mathrm{\nabla }\cdot \hat{\mathit{\sigma }}(𝐱,\omega )=-{\omega }^2\hat{\rho }(𝐱,\omega )\widehat{𝐮}(𝐱,\omega ).$

Substituting equation (5) into equation (8) we get

$\text{(9)}\mathrm{\nabla }\cdot (𝖢:\widehat{\mathit{ϵ}})+{\omega }^2\hat{\rho }\widehat{𝐮}=0.$

Also, taking the Fourier transform of equation (3), we have

$\text{(10)}\widehat{\mathit{ϵ}}=\frac{1}{2}[\mathrm{\nabla }\widehat{𝐮}+(\mathrm{\nabla }\widehat{𝐮}{)}^T].$

Since $\hat{\mathit{\sigma }}$ and $\widehat{\mathit{ϵ}}$ are symmetric, we must have

$C_{ijkl}=C_{jikl}=C_{ijlk}.$

Because of this symmetry, we can replace $\widehat{\mathit{ϵ}}$ by $\mathrm{\nabla }\widehat{𝐮}$ in equation (9) to get

$\mathrm{\nabla }\cdot (𝖢:\mathrm{\nabla }\widehat{𝐮}\mathit{)}+{\omega }^2\hat{\rho }\widehat{𝐮}=0.$

Dropping the hats, we then get the wave equation for elastodynamics

$\text{(11)}\mathrm{\nabla }\cdot (𝖢:\mathrm{\nabla }𝐮\mathit{)}+{\omega }^2\rho 𝐮=0.$

Antiplane shear

Let us now consider the case of antiplane shear. Assume that $𝖢$ is isotropic, i.e.,

$𝖢:\mathrm{\nabla }𝐮=\mu [\mathrm{\nabla }𝐮+(\mathrm{\nabla }𝐮{)}^T]+\lambda \text{tr}(\mathrm{\nabla }𝐮)1$

where $\mu$ is the shear modulus and $\lambda$ is the Lame modulus.

Let us assume that $\mu$ and $\lambda$ are independent of $x_1$, i.e.,

$\mu \equiv \mu (x_2,x_3)\text{and}\lambda \equiv \lambda (x_2,x_3).$

Let us look for a solution with $u_2=u_3=0$ and $u_1$ independent of $x_1$, i.e., $u_1\equiv u_1(x_2,x_3)$. This is an out of plane mode of deformation.

Then, noting that $\text{tr}(\mathrm{\nabla }𝐮)=\mathrm{\nabla }\cdot 𝐮$, we have

$\text{tr}(\mathrm{\nabla }𝐮)=\frac{\partial u_i}{\partial x_i}=0.$

Therefore,

$[𝖢:\mathrm{\nabla }𝐮{]}_{ij}=\mu [\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}]$

or

$𝖢:\mathrm{\nabla }𝐮=\mu \left[\begin{array}{ccc}0& \frac{\partial u_1}{\partial x_2}& \frac{\partial u_1}{\partial x_3}\\ \\ \frac{\partial u_1}{\partial x_2}& 0& 0\\ \\ \frac{\partial u_1}{\partial x_3}& 0& 0\end{array}\right].$

Therefore

$\mathrm{\nabla }\cdot (𝖢:\mathrm{\nabla }𝐮\mathit{)}=\left[\begin{array}{c}\frac{\partial }{\partial x_2}\left(\mu \frac{\partial u_1}{\partial x_2}\right)+\frac{\partial }{\partial x_3}\left(\mu \frac{\partial u_1}{\partial x_3}\right)\\ 0\\ 0\end{array}\right].$

Plugging into the wave equation (11) we get

$\frac{\partial }{\partial x_2}\left(\mu \frac{\partial u_1}{\partial x_2}\right)+\frac{\partial }{\partial x_3}\left(\mu \frac{\partial u_1}{\partial x_3}\right)+{\omega }^2\rho u_1=0$

or (using the two-dimensional gradient operator $\bar{\mathrm{\nabla }}$)

$\text{(12)}\bar{\mathrm{\nabla }}\cdot (\mu \overline{\mathrm{\nabla }{u}_1)}+{\omega }^2\rho u_1=0.$

TM and TE modes in electromagnetism

Let us now consider the TM (transverse magnetic field) and TE (transverse electric field) modes in electromagnetism and look for parallels with antiplane shear in elastodynamics.

Recall the Maxwell equations (with hats dropped)

$\text{(13)}\mathrm{\nabla }\times 𝐄=i\omega \mu 𝐇;\mathrm{\nabla }\times 𝐇=-i\omega ϵ𝐄.$

Assume that $\mu$ and $ϵ$ are scalars which are independent of $x_1$, i.e., $\mu \equiv \mu (x_2,x_3)$ and $ϵ\equiv ϵ(x_2,x_3)$.

For the TE case, we look for solutions with $E_2=E_3=0$ and $E_1$ independent of $x_1$, i.e., $E_1\equiv E_1(x_2,x_3)$.

Then,

$\mathrm{\nabla }\times 𝐄=[0,\frac{\partial E_1}{\partial x_3},-\frac{\partial E_1}{\partial x_2}].$

This implies that

$𝐇=[0,\frac{1}{i\omega \mu }\frac{\partial E_1}{\partial x_3},-\frac{1}{i\omega \mu }\frac{\partial E_1}{\partial x_2}].$

Therefore,

$\mathrm{\nabla }\times 𝐇=[-\frac{\partial }{\partial x_3}\left(\frac{1}{i\omega \mu }\frac{\partial E_1}{\partial x_3}\right)-\frac{\partial }{\partial x_2}\left(\frac{1}{i\omega \mu }\frac{\partial E_1}{\partial x_2}\right),0,0].$

or,

$\mathrm{\nabla }\times 𝐇=\frac{i}{\omega }[\frac{\partial }{\partial x_2}\left(\frac{1}{\mu }\frac{\partial E_1}{\partial x_2}\right)+\frac{\partial }{\partial x_3}\left(\frac{1}{\mu }\frac{\partial E_1}{\partial x_3}\right),0,0]=\frac{i}{\omega }[\bar{\mathrm{\nabla }}\cdot \left(\frac{1}{\mu (x_2,x_3)}\bar{\mathrm{\nabla }}{E}_1\right),0,0].$

Plugging into equation (13) we get the TE equation

$\bar{\mathrm{\nabla }}\cdot \left(\frac{1}{\mu (x_2,x_3)}\bar{\mathrm{\nabla }}{E}_1\right)+{\omega }^2ϵ(x_2,x_3)E_1=0.$

This equation has the same form as equation (12).

More generally, if

$\mathit{\mu }=\mathit{\mu }(x_2,x_3)=\left[\begin{array}{ccc}{\mu }_{11}& 0& 0\\ 0& {\mu }_{22}& {\mu }_{23}\\ 0& {\mu }_{23}& {\mu }_{33}\end{array}\right]=\left[\begin{array}{cc}{\mu }_{11}& \left[\mathrm{𝟢}\right]\\ \left[\mathrm{𝟢}\right]& \left[𝖬\right]\end{array}\right]\text{where}\left[𝖬\right]=\left[\begin{array}{cc}{\mu }_{22}& {\mu }_{23}\\ {\mu }_{23}& {\mu }_{33}\end{array}\right]\equiv 𝑴$

and

$\mathit{ϵ}=\mathit{ϵ}(x_2,x_3)=\left[\begin{array}{cc}{ϵ}_{11}& \left[\mathrm{𝟢}\right]\\ \left[\mathrm{𝟢}\right]& \left[𝖭\right]\end{array}\right]\text{where}\left[𝖭\right]=\left[\begin{array}{cc}{ϵ}_{22}& {ϵ}_{23}\\ {ϵ}_{23}& {ϵ}_{33}\end{array}\right]\equiv 𝑵$

we get the TE equation

$\bar{\mathrm{\nabla }}\cdot [({𝑹}_{\perp }\cdot {𝑴}^{-1}\cdot {𝑹}_{\perp }^T)\cdot \bar{\mathrm{\nabla }}{E}_1]+{\omega }^2{ϵ}_{11}{E}_1\mathrm{𝟏}=\mathrm{𝟎}\text{where}{𝑹}_{\perp }\equiv {\left[𝖱\right]}_{\perp }=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right].$

Similarly, there is a TM equation with $H_2=H_3=0$ of the form

$\bar{\mathrm{\nabla }}\cdot [({𝑹}_{\perp }\cdot {𝑵}^{-1}\cdot {𝑹}_{\perp }^T)\cdot \bar{\mathrm{\nabla }}{H}_1]+{\omega }^2{\mu }_{11}{H}_1\mathrm{𝟏}=\mathrm{𝟎}$

which for the isotropic case reduces to

$\bar{\mathrm{\nabla }}\cdot \left(\frac{1}{ϵ(x_2,x_3)}\bar{\mathrm{\nabla }}{H}_1\right)+{\omega }^2\mu (x_2,x_3)H_1=0.$

The general solution independent of $x_1$ is a superposition of the TE and TM solutions. This can be seen by observing that the Maxwell equations decouple under these conditions and a general solution can be written as

$(E_1,E_2,E_3)=(E_1,0,0)+(0,E_2,E_3)$

where the first term represents the TE solution. We can show that the second term represents the TM solution by observing that

$\mathrm{\nabla }\times (0,E_2,E_3)=[\frac{\partial E_3}{\partial x_2}-\frac{\partial E_2}{\partial x_3},0,0]$

implying that $H_2=H_3=0$ which is the TM solution.

A resonant structure

Consider the periodic geometry shown in Figure 2. The matrix material has a high value of shear modulus ($\mu$) while the split-ring shaped region has a low shear modulus or is a void. The material inside the ring has the same shear modulus as the matrix material and is connected to the matrix by a thin ligament. The system is subjected to a displacement $u_1$ in the $x_1$ direction (parallel to the axis of each cylindrical split ring).

File:SplitCylArray.jpg Figure 2. A periodic geometry containing split hollow cylinders of soft material in a matrix of stiff material. The $x_1$ direction is parallel to the axis of each cylinder.

Clearly, each periodic component of the system behaves like a mass attached to a spring. This is a resonant structure and the effective density ${\rho }_{11}^{*}(\omega )$ can be negative. A detailed treatment of the problem can be found in Movchan04. Note that the governing equation for this problem is

$\bar{\mathrm{\nabla }}\cdot (\mu \bar{\mathrm{\nabla }}{u}_1)+{\omega }^2\rho u_1=0.$

Let us compare this problem with the TM case where $H_1$ is the out of plane magnetic induction. The governing equation now is

$\bar{\mathrm{\nabla }}\cdot (\frac{1}{ϵ}\bar{\mathrm{\nabla }}{H}_1)+{\omega }^2\mu H_1=0.$

If the value of $1/ϵ$ in the region of the void (ring) is small and hence $ϵ$ is large (which implies that the conductivity $\sigma$ is large), analogy with the equation of elastodynamics implies that the effective permeability ${\mu }_{11}^{*}$ can be negative for this material.

References

• R. J. Atkin and N. Fox. An introduction to the theory of elasticity. Longman, New York, 1980.
• A. B. Movchan and S. Guenneau. Split-ring resonators and localized modes. Physical Review B, 70:125116, 2004.