The energy density for an electromagnetic field is given by
$$ u=\frac 12 \epsilon_0 |\vec E|^2+\frac{1}{2\mu_0}|\vec B|^2 $$
And the energy $U$ is calculated by integrating the energy density over all space
$$ U=\int\text dV \,u $$
We can use the electric and magnetic fields which describe plane waves to compute the energy density
$$ \begin{aligned} u&=\frac 12 \epsilon_0 E^2_0\left [ \mathrm{Re}\left ( e^{i(\vec k \cdot r-\omega t+\theta)} \right ) \right ]^2+\frac{1}{2\mu_0 c^2}E^2_0\left [ \mathrm{Re}\left ( e^{i(\vec k \cdot r-\omega t+\theta)} \right ) \right ]^2 \\ &=\epsilon_0 E^2_0 \left [ \mathrm{Re}\left ( e^{i(\vec k \cdot r-\omega t+\theta)} \right ) \right ]^2
\end{aligned} $$
Notice that the contributions to $u$ from the electric and magnetic fields are identical.
We can read off that the maximum energy density is
$$ u_\text{max}=\epsilon_0 E^2_0 $$
and we can calculate that the average energy density is
$$ \left < u \right > =\frac 12 \epsilon_0 E^2_0 $$
as $\left < \cos^2 (\vec k \cdot \vec r-\omega t+\theta) \right >=\frac 12$
💼 Case: consider the volume covered by a wave in the time-interval $\Delta t$ through the surface are $\Delta A$
The average energy density enclosed by this volume is given by
$$ \left < \text{enclosed energy density} \right > = \left < u \right > c\Delta t\Delta A=\frac 12 \epsilon_0 E^2_0 c\Delta t \Delta A=\Delta U_\text{tot} $$
Hence the average flux through a surface is given by
$$ \mathcal F=\frac{\Delta U_\text{tot}}{\Delta t\Delta A}=\frac{E^2_0}{2Z_0} $$
where we have defined the impedance of free space
$$ Z_0\equiv \sqrt{\frac{\mu_0}{\epsilon_0}}=377\, \Omega $$
We now calculate the rate of the energy density
$$ \dot u =\frac 12 \epsilon_0 \frac{\partial}{\partial t} ( \vec E\cdot \vec E)+\frac{1}{2\mu_0}\frac{\partial}{\partial t}(\vec B\cdot \vec B)=\epsilon_0 \vec E\cdot \dot{\vec E} +\frac 1 {\mu_0} \vec B\cdot \dot{\vec B} $$
Substituting with Maxwell equations in free space $\dot{\vec B}=-\vec \nabla \times \vec E$ , $\dot{\vec E}=c^2 \vec \nabla \times \vec B$ we obtain
$$ \dot u=\frac{1}{\mu_0} (\vec E\cdot \vec \nabla \times \vec B-\vec B\cdot \vec \nabla \times \vec E)=-\vec \nabla \cdot N $$
where we have defined
$$ \vec N=\frac{1}{\mu_0} \vec E\times \vec B $$
which is the Poynting vector
It is an energy flux, the units of the Poynting vector are
$$ \rm [N]=\frac{V\,m^{-1} T}{Wb \, A^{-1} \, m^{-1}}=A \, V\, m^{-2}=W\, m^{-2} $$
Watts per square meter
The rate of change of energy is
$$ \dot U=-\int \text dV \,\vec \nabla \cdot \vec N =-\int \vec N \cdot \text d\vec S $$
where the second equality arises from the divergence theorem
💎 Conclusions:
We identify $\vec N$ as the energy flux through a surface
In materials it is often useful to define the Poynting vector as
$$ \vec N =\vec E\times \vec H $$
which differs to the earlier definition if $\mu_r\ne 1$. The difference arises due to the presence of currents ( and energy loss to driving these currents). In the former definition the energy expended driving all currents is considered. In latter $\vec E\times \vec H$ the energy expended on $free$ currents alone is balanced by the flux
Using $\vec N=\frac{1}{\mu_0} \vec E\times \vec B$ we can insert a plane wave solution $\vec E=\vec E_0 e^{i(\vec k\cdot \vec r-\omega t)}$ and $\vec B=\vec B_0 e^{i(\vec k\cdot \vec r-\omega t)}$
$$ \begin{aligned} \vec N&=\frac{1}{\mu_0}\vec E\times \vec B =\frac{\vec E_0 \times (\vec k \times \vec E_0)}{\omega \mu_0}\left [ \mathrm{Re}(e^{i(\vec k \cdot \vec r-\omega t+\phi)}) \right ]^2 \\ &=\frac{E^2_0}{Z_0}\hat k \cos^2(\vec k \cdot \vec r-\omega t +\theta)
\end{aligned} $$
where $Z_0$ is the impedance of free space.
🗒️ Note: that $\hat k$ is in the same direction as $\vec N$ so that the waves energy flux goes in the same direction as the wavevector (which is not the case for the electric or magnetic fields).
The maximum value that the Poynting vector takes is
$$ N_\text{max}=\frac{E_0^2}{Z_0} $$
🗒️ Note: that the units of $N_\text{max}$ are $\rm W\, m^{-2}$.
The irradiance ( or intensity) is defined to be the average of the Poynting vector
$$ \left < N \right > =\frac{E^2_0}{2 Z_0} $$
The momentum of a photon is $p=U/c$ and hence the momentum flux is $\left < N \right >/c$ [$\rm Nm^{-2}$], and this is the radiation pressure, $P_r$.
At a boundary the radiation creates a force $F_r$ and the magnitude of this force depends on whether or not the radiation is absorbed or reflected by the boundary
If all radiation is absorbed then
$$ F_r=\frac{A\left < N \right >}{c} $$
If all radiation is reflected then
$$ F_r=\frac{2A\left < N \right >}{c} $$