The D'Alemnertian of the Electric Field Inside a Dielectric With No Free Charges

To find the D'Alembertian of the electric field  
\[\mathbf{E}\]
  in a dieletric with no free charges, start with the Maxwell equations for a dielectric.
\[\mathbf{\nabla} \cdot \mathbf{E} =\frac{\mathbf{\nabla} \cdot \mathbf{p}}{\epsilon_0}\]
 (1)
\[\mathbf{\nabla} \times \mathbf{E} =- \frac{\partial \mathbf{H}}{\partial t}\]
 (2)
\[c^2 \mathbf{\nabla} \times \mathbf{H} = \frac{\partial}{\partial t} ( \frac{1}{\epsilon_o}\mathbf{p} + \mathbf{E}\]
 (3)
\[\mathbf{\nabla} \cdot \mathbf{H}=0\]
 (4)
where  
\[\mathbf{p}\]
  is the polarization vector. The D'Alembertian operator is  
\[\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} -\frac{1}{c^2} \frac{\partial^2}{\partial t^2} = \nabla^2 - \frac{1}{c^2} \frac{\partial^2}{\partial t^2}\]
&
Start by taking the curl of (2) to give
\[\mathbf{\nabla} \times (\mathbf{\nabla} \times \mathbf{E}) =- \mathbf{\nabla} \times ( \frac{\partial \mathbf{H}}{\partial t})\]

Use the vector identity  
\[\mathbf{\nabla} \times (\mathbf{\nabla} \times \mathbf{E}) =\mathbf{\nabla} (\mathbf{\nabla} \cdot \mathbf{E}) - \nabla^2 \mathbf{E}\]
  to rewrite this as
\[\mathbf{\nabla} (\mathbf{\nabla} \cdot \mathbf{E}) - \nabla^2 \mathbf{E} =- \mathbf{\nabla} \times ( \frac{\partial \mathbf{H}}{\partial t})\]

Substitute (3) into this to obtain
\[\begin{equation} \begin{aligned} \mathbf{\nabla} (\mathbf{\nabla} \cdot \mathbf{E}) - \nabla^2 \mathbf{E}& = - \frac{\partial^2 }{\partial t^2 } ( \frac{1}{c^2 \epsilon_o} \mathbf{p} + \frac{\mathbf{E}}{c^2}) \\ &= - \frac{1}{c^2 \epsilon_o} \frac{\partial^2 \mathbf{p}}{\partial t^2} - \frac{\partial^2 \mathbf{E}}{c^2 \partial t^2} \end{aligned} \end{equation}\]

From (1) we obtain
\[\mathbf{\nabla} (\mathbf{\nabla} \cdot \mathbf{E}) =\mathbf{\nabla} (\frac{\mathbf{\nabla} \cdot \mathbf{p}}{\epsilon_0})=\frac{1}{\epsilon_0} \mathbf{\nabla} ({\mathbf{\nabla} \cdot \mathbf{p}})\]
 (1)
Substituting this equation into the previous one and simplifying gives>
\[\nabla^2 \mathbf{E} - \frac{1}{c^2} \frac{\partial^2 \mathbf{E}}{\partial t^2} =- \frac{1}{\epsilon_0} \mathbf{\nabla} ( \mathbf{\nabla} \cdot \mathbf{E}) + \frac{1}{c^2 \epsilon_0} \frac{\partial^2 \mathbf{p}}{\partial t^2 } \]