真空の場合のマクスウェル方程式
簡単のために真空（$\rho = 0, \ \boldsymbol{J} = \boldsymbol{0}$）の場合で説明する。マクスウェル方程式は，$ \boldsymbol{D} = \varepsilon_0 \boldsymbol{E}, \ \boldsymbol{H} = \varepsilon_0 c^2 \boldsymbol{B}$ を使って $\boldsymbol{E}$ と $\boldsymbol{B}$ であらわすと，
\begin{eqnarray}\nabla\cdot \boldsymbol{E} &=& 0 \tag{1}\\
\nabla\cdot\boldsymbol{B} &=& 0 \tag{2}\\
\frac{\partial \boldsymbol{B}}{\partial t} &=& – \nabla\times\boldsymbol{E} \tag{3}\\
\frac{\partial \boldsymbol{E}}{\partial t} &=&c^2 \nabla\times\boldsymbol{B} \tag{4}
\end{eqnarray}
$(3)$ 式と $(4)$ 式から
\begin{eqnarray}
\frac{\partial}{\partial t} \left\{\frac{1}{2} \varepsilon_0 \boldsymbol{E}\cdot\boldsymbol{E} + \frac{1}{2} \varepsilon_0 c^2 \boldsymbol{B}\cdot\boldsymbol{B} \right\}
&=&
\varepsilon_0 \boldsymbol{E}\cdot\frac{\partial \boldsymbol{E}}{\partial t} + \varepsilon_0 c^2 \boldsymbol{B}\cdot\frac{\partial \boldsymbol{B}}{\partial t} \\
&=& \varepsilon_0 \boldsymbol{E}\cdot\left(c^2 \nabla\times\boldsymbol{B} \right) + \varepsilon_0 c^2 \boldsymbol{B}\cdot\left( – \nabla\times\boldsymbol{E} \right) \\
&=& – \varepsilon_0 c^2 \left\{ \left( \nabla\times\boldsymbol{E} \right)\cdot\boldsymbol{B} – \boldsymbol{E} \cdot \left( \nabla\times\boldsymbol{B} \right) \right\} \\
&=& – \varepsilon_0 c^2 \nabla\cdot \left(\boldsymbol{E}\times\boldsymbol{B} \right)
\end{eqnarray}
（最後の行になるのはこのへんを参照。）
連続の式：エネルギー保存則
あらためて
\begin{eqnarray}
u &\equiv& \frac{1}{2} \varepsilon_0 \boldsymbol{E}\cdot\boldsymbol{E} + \frac{1}{2} \varepsilon_0 c^2 \boldsymbol{B}\cdot\boldsymbol{B} \\
&=& \frac{1}{2}\boldsymbol{E}\cdot\boldsymbol{D} + \frac{1}{2} \boldsymbol{H}\cdot\boldsymbol{B}\\
\boldsymbol{S} &\equiv& \varepsilon_0 c^2\boldsymbol{E}\times\boldsymbol{B} \\
&=& \boldsymbol{E}\times\boldsymbol{H}
\end{eqnarray}
と定義すれば，
$$\frac{\partial u}{\partial t} + \nabla\cdot\boldsymbol{S} = 0$$
という連続の式が出てくる。別ページ「参考：電荷の保存則」で説明したように，連続の式というのは保存則の微分形である。
これはすなわち，$u$ が電磁場のエネルギー密度，$\boldsymbol{S}$ が電磁場のエネルギー流束をあらわすことを意味する。電磁場自体がもつエネルギーの流れをあらわすベクトル $\boldsymbol{S}$ はポインティングベクトルと呼ばれる。
いきなり，$\displaystyle u = \frac{1}{2}\boldsymbol{E}\cdot\boldsymbol{D} + \frac{1}{2} \boldsymbol{H}\cdot\boldsymbol{B}$ が電磁場のエネルギー密度だと言われても唐突で困るだろうから，以下ではまず，静電場の場合に，$\displaystyle u = \frac{1}{2}\boldsymbol{E}\cdot\boldsymbol{D}$ が確かに静電場のエネルギー密度になっていることを確かめてみる。
静電場のエネルギー
2個の荷電粒子の静電ポテンシャルエネルギー
位置 $\boldsymbol{r}_1$ に電荷 $q_1$，位置 $\boldsymbol{r}_2$ に電荷 $q_2$。この1対の荷電粒子の静電ポテンシャルエネルギー $U$ は，電荷 $q_2$ が位置 $\boldsymbol{r}_1$ につくる静電ポテンシャルを $\phi_2(\boldsymbol{r}_1)$ とすると，
$$U = U_{12} = q_1 \phi_2(\boldsymbol{r}_1)$$
あるいは，
電荷 $q_1$ が位置 $\boldsymbol{r}_2$ につくる静電ポテンシャルを $\phi_1(\boldsymbol{r}_2)$ とすると，
$$U = U_{21} = q_2 \phi_1(\boldsymbol{r}_2)$$
2個で1対のペア1組に対してポテンシャルエネルギーがある。1個1個の荷電粒子について機械的に計算する場合は，ダブらないように $2$ で割る。
$$U = \frac{1}{2} \left(U_{12} + U_{21}\right)$$
$U_{11}$ とか $U_{22}$ とかは無し。
3個の荷電粒子の系のポテンシャルエネルギー
\begin{eqnarray}
U &=& U_{12} + U_{13} + U_{23} \\
&=& \frac{1}{2} \left\{U_{12} + U_{13} + U_{21} + U_{23} + U_{31} + U_{32} \right\} \\
&=& \frac{1}{2} \sum_{i=1}^3 q_i \phi(\boldsymbol{r}_i) \\
\phi(\boldsymbol{r}_i) &\equiv& \sum_{j\neq i} \phi_j (\boldsymbol{r}_i)
\end{eqnarray}
$n$ 個の荷電粒子の静電ポテンシャルエネルギー
同様にして，荷電粒子が $n$ 個の場合は
\begin{eqnarray}
U
&=& \frac{1}{2} \sum_{i=1}^n q_i \phi(\boldsymbol{r}_i)
\end{eqnarray}
連続的な電荷密度の場合の静電ポテンシャルエネルギー
連続的な電荷分布 \(\rho(\boldsymbol{r})\) の場合は，以下のような置き換えをする。
\begin{eqnarray}
\boldsymbol{r}_i &\rightarrow& \boldsymbol{r} \\
q_i &\rightarrow& \rho(\boldsymbol{r}_i ) dV_i\\
\sum_{i} dV_i&\rightarrow& \iiint \, dV
\end{eqnarray}
また，$\nabla\cdot\boldsymbol{D}=\rho, \ \boldsymbol{E} = – \nabla\phi$ も使って
\begin{eqnarray}
U &=& \frac{1}{2} \sum_{i=1}^n q_i \phi(\boldsymbol{r}_i) \\
&\rightarrow& \frac{1}{2} \iiint_V \rho \,\phi\, dV
\end{eqnarray}
つまり，静電場の場合の静電ポテンシャルのエネルギー密度 $u_{\rm e.s.}$ は
$$u_{\rm e.s.} = \frac{1}{2}\rho \,\phi$$
と書けることがわかった。
静電場のエネルギー密度
電場のエネルギーが
$$ U = \frac{1}{2} \iiint_V \boldsymbol{E}\cdot\boldsymbol{D}\, dV$$
と書けるとすると，静電場の場合には
\begin{eqnarray}
U &=& \frac{1}{2} \iiint_V \boldsymbol{E}\cdot\boldsymbol{D}\, dV \\
&=& – \frac{1}{2} \iiint_V \left(\nabla \phi\right)\cdot\boldsymbol{D}\, dV \\
&=& – \frac{1}{2} \iiint_V \left\{ \nabla\cdot \left( \phi \boldsymbol{D}\right) -\left(\nabla\cdot\boldsymbol{D}\right) \phi \right\} \, dV \\
&=& – \frac{1}{2} \iint_S \left( \phi \boldsymbol{D}\right)\cdot\boldsymbol{n}\, dS
+ \frac{1}{2} \iiint_V \rho\,\phi\, dV \\
&\rightarrow & \frac{1}{2} \iiint_V \rho\,\phi\, dV
\end{eqnarray}
ここで，ガウスの定理によって変形した表面積分の項 $\displaystyle \iint_S \left(\phi \boldsymbol{D} \right)\cdot\boldsymbol{n} \,dS$ は十分大きい領域 $S$ を考えればゼロになることを使った。
よって，電場のエネルギー密度 $\displaystyle u = \frac{1}{2} \boldsymbol{E}\cdot \boldsymbol{D}$ は，静電場の場合には静電ポテンシャルのエネルギー密度 $u_{\rm e.s.}$ に帰着することが示された。
静磁場のエネルギー密度
では，静磁場の場合に，磁場のエネルギー密度はどのように書けるだろうか。
磁場のエネルギーが
$$ U = \frac{1}{2} \iiint_V \boldsymbol{H}\cdot\boldsymbol{B}\, dV$$
と書けるとすると，静磁場の場合には
\begin{eqnarray}
U &=& \frac{1}{2} \iiint_V \boldsymbol{H}\cdot\boldsymbol{B}\, dV \\
&=& \frac{1}{2} \iiint_V \boldsymbol{H}\cdot\left( \nabla\times\boldsymbol{A}\right)\, dV \\
&=& \frac{1}{2} \iiint_V \left\{ \left( \nabla\times\boldsymbol{H} \right)\cdot\boldsymbol{A}
-\nabla\cdot\left(\boldsymbol{H}\times\boldsymbol{A} \right)\right\} \, dV \\
&=& \frac{1}{2} \iiint_V \boldsymbol{J}\cdot\boldsymbol{A} \, dV
-\frac{1}{2} \iint_S \left(\boldsymbol{H}\times\boldsymbol{A} \right)\cdot \boldsymbol{n}\, dS \\
&\rightarrow& \frac{1}{2} \iiint_V \boldsymbol{J}\cdot\boldsymbol{A} \, dV
\end{eqnarray}
ここで，ガウスの定理によって変形した表面積分の項 $\displaystyle \iint_S \left(\boldsymbol{H}\times\boldsymbol{A} \right)\cdot \boldsymbol{n} \,dS$ は十分大きい領域 $S$ を考えればゼロになることを使った。
よって，静磁場の場合のエネルギー密度 $u_{\rm m.s.}$ は
$$u_{\rm m.s.} = \frac{1}{2} \boldsymbol{J}\cdot\boldsymbol{A}$$
と書けることが示された。静電場のエネルギー密度 $u_{\rm e.s.}$ が電荷密度 $\rho$ とスカラーポテンシャル $\phi$ の積で書かれることに対応して，静磁場のエネルギー密度 $u_{\rm m.s.}$ は電流密度 $\boldsymbol{J}$ とベクトルポテンシャル $\boldsymbol{A}$ の内積で書かれていることがわかって，なんだかとても清々しい。
補足：ベクトルの微分公式
以下のようなベクトルの微分の公式を使っている。（右辺から左辺を引いてゼロになることを直接示せる。こことかこことかを参照。）
$$\nabla\cdot\left( \boldsymbol{E}\times\boldsymbol{B}\right)
= \left(\nabla\times\boldsymbol{E} \right)\cdot\boldsymbol{B}
– \boldsymbol{E}\cdot\left(\nabla\times\boldsymbol{B} \right)$$
$$\nabla\cdot\left( \phi\,\boldsymbol{D}\right) = \left(\nabla\phi\right)\cdot\boldsymbol{D} +\left(\nabla\cdot\boldsymbol{D}\right)\,\phi$$