In the frequency domain

In the frequency domain with $$\boldsymbol{\mathcal{E}} = \textrm{Re} \BE e^{j \omega t}, \boldsymbol{\mathcal{H}} = \textrm{Re} \BH e^{j \omega t}$$. Using the electric field dot product as an example, note that we can write

\label{eqn:energyMomentumWithMagneticSources:480}
\boldsymbol{\mathcal{E}} = \inv{2} \lr{ \BE e^{j \omega t} + \BE^\conj e^{-j \omega t} },

so

\label{eqn:energyMomentumWithMagneticSources:500}
\begin{aligned}
\boldsymbol{\mathcal{E}}^2
&=
\inv{2} \lr{ \BE e^{j \omega t} + \BE^\conj e^{-j \omega t} }
\cdot
\inv{2} \lr{ \BE e^{j \omega t} + \BE^\conj e^{-j \omega t} } \\
&=
\inv{4} \lr{
\BE^2 e^{2 j \omega t}
+ \BE \cdot \BE^\conj + \BE^\conj \cdot \BE
+\lr{\BE^\conj}^2 e^{-2 j \omega t}
} \\
&=
\inv{2} \textrm{Re}
\lr{
\BE \cdot \BE^\conj
+
\BE^2 e^{2 j \omega t}
}.
\end{aligned}

Similarly, for the cross product

\label{eqn:energyMomentumWithMagneticSources:540}
\begin{aligned}
\boldsymbol{\mathcal{E}} \cross \boldsymbol{\mathcal{H}}
&=
\inv{4}
\lr{
\BE \cross \BH e^{2 j \omega t}
+ \BE \cross \BH^\conj + \BE^\conj \cross \BH
+ \lr{ \BE^\conj \cross \BH^\conj } e^{-2 j \omega t}
} \\
&=
\inv{2}
\textrm{Re}
\lr{
\BE \cross \BH^\conj
+
\BE \cross \BH e^{2 j \omega t}
}.
\end{aligned}

Given phasor representations of the sources $$\boldsymbol{\mathcal{M}} = \BM e^{j \omega t}, \boldsymbol{\mathcal{J}} = \BJ e^{j \omega t}$$, \ref{eqn:energyMomentumWithMagneticSources:40} can be recast into (a messy) phasor form

\label{eqn:energyMomentumWithMagneticSources:560}
\begin{aligned}
\inv{2} &\textrm{Re} \inv{2} \PD{t}{} \lr{
\epsilon_0 \BE \cdot \BE^\conj
+ \mu_0 \BH \cdot \BH^\conj
+ \epsilon_0 \BE^2 e^{ 2 j \omega t}
+ \mu_0 \BH^2 e^{ 2 j \omega t}
} \\
&+
\BE \cross \BH^\conj
+\BE \cross \BH e^{ 2 j \omega t}
} \\
&=
\inv{2} \textrm{Re}
\lr{
– \BH \cdot \BM^\conj
– \BE \cdot \BJ^\conj
– \BH \cdot \BM e^{2 j \omega t}
– \BE \cdot \BJ e^{2 j \omega t}
}.
\end{aligned}

In particular, when averaged over one period, the oscillatory terms vanish. The time averaged equivalent of the Poynting theorem is thus

\label{eqn:energyMomentumWithMagneticSources:580}
0 =
{\left[
\textrm{Re}
\lr{
\inv{2} \PD{t}{} \lr{
\epsilon_0 \BE \cdot \BE^\conj
+ \mu_0 \BH \cdot \BH^\conj
}
+
\BE \cross \BH^\conj
}
+
\BH \cdot \BM^\conj
+
\BE \cdot \BJ^\conj
}
\right]
}_{\textrm{av}}.

Comparison to the reciprocity theorem result

The reciprocity theorem had a striking similarity to the Poynting theorem above, which isn’t suprising since both were derived by calculating the divergence of a Poynting like quantity.

Here’s a repetition of the reciprocity divergence calculation without the single frequency (phasor) assumption

\label{eqn:energyMomentumWithMagneticSources:600}
\begin{aligned}
\boldsymbol{\mathcal{E}}^{(a)} \cross \boldsymbol{\mathcal{H}}^{(b)}
-\boldsymbol{\mathcal{E}}^{(b)} \cross \boldsymbol{\mathcal{H}}^{(a)}
} \\
&=
\boldsymbol{\mathcal{H}}^{(b)} \cdot \lr{ \spacegrad \cross \boldsymbol{\mathcal{E}}^{(a)} } -\boldsymbol{\mathcal{E}}^{(a)} \cdot \lr{ \spacegrad \cross \boldsymbol{\mathcal{H}}^{(b)} } \\
-\boldsymbol{\mathcal{H}}^{(a)} \cdot \lr{ \spacegrad \cross \boldsymbol{\mathcal{E}}^{(b)} } +\boldsymbol{\mathcal{E}}^{(b)} \cdot \lr{ \spacegrad \cross \boldsymbol{\mathcal{H}}^{(a)} } \\
&=
-\boldsymbol{\mathcal{H}}^{(b)} \cdot \lr{ \mu_0 \partial_t \boldsymbol{\mathcal{H}}^{(a)} + \boldsymbol{\mathcal{M}}^{(a)} }
-\boldsymbol{\mathcal{E}}^{(a)} \cdot \lr{ \boldsymbol{\mathcal{J}}^{(b)} + \epsilon_0 \partial_t \boldsymbol{\mathcal{E}}^{(b)} } \\
+\boldsymbol{\mathcal{H}}^{(a)} \cdot \lr{ \mu_0 \partial_t \boldsymbol{\mathcal{H}}^{(b)} + \boldsymbol{\mathcal{M}}^{(b)} }
+\boldsymbol{\mathcal{E}}^{(b)} \cdot \lr{ \boldsymbol{\mathcal{J}}^{(a)} + \epsilon_0 \partial_t \boldsymbol{\mathcal{E}}^{(a)} } \\
&=
\epsilon_0
\lr{
\boldsymbol{\mathcal{E}}^{(b)} \cdot \partial_t \boldsymbol{\mathcal{E}}^{(a)}
-\boldsymbol{\mathcal{E}}^{(a)} \cdot \partial_t \boldsymbol{\mathcal{E}}^{(b)}
}
+
\mu_0
\lr{
\boldsymbol{\mathcal{H}}^{(a)} \cdot \partial_t \boldsymbol{\mathcal{H}}^{(b)}
-\boldsymbol{\mathcal{H}}^{(b)} \cdot \partial_t \boldsymbol{\mathcal{H}}^{(a)}
} \\
&+\boldsymbol{\mathcal{H}}^{(a)} \cdot \boldsymbol{\mathcal{M}}^{(b)}
-\boldsymbol{\mathcal{H}}^{(b)} \cdot \boldsymbol{\mathcal{M}}^{(a)}
+\boldsymbol{\mathcal{E}}^{(b)} \cdot \boldsymbol{\mathcal{J}}^{(a)}
-\boldsymbol{\mathcal{E}}^{(a)} \cdot \boldsymbol{\mathcal{J}}^{(b)}
\end{aligned}

What do these time derivative terms look like in the frequency domain?

\label{eqn:energyMomentumWithMagneticSources:620}
\begin{aligned}
\boldsymbol{\mathcal{E}}^{(b)} \cdot \partial_t \boldsymbol{\mathcal{E}}^{(a)}
&=
\inv{4}
\lr{
\BE^{(b)} e^{j \omega t}
+
{\BE^{(b)}}^\conj e^{-j \omega t}
}
\cdot
\partial_t
\lr{
\BE^{(a)} e^{j \omega t}
+
{\BE^{(a)}}^\conj e^{-j \omega t}
} \\
&=
\frac{j \omega}{4}
\lr{
\BE^{(b)} e^{j \omega t}
+
{\BE^{(b)}}^\conj e^{-j \omega t}
}
\cdot
\lr{
\BE^{(a)} e^{j \omega t}

{\BE^{(a)}}^\conj e^{-j \omega t}
} \\
&=
\frac{\omega}{4}
\lr{
j \BE^{(a)} \cdot { \BE^{(b)} }^\conj
-j \BE^{(b)} \cdot { \BE^{(a)} }^\conj
+j \BE^{(a)} \cdot \BE^{(b)} e^{ 2 j \omega t }
-j { \BE^{(a)}}^\conj \cdot { \BE^{(b)} }^\conj e^{ -2 j \omega t }
} \\
&=
\inv{2} \textrm{Re}
\lr{
j \omega \BE^{(a)} \cdot { \BE^{(b)} }^\conj
+ j \omega \BE^{(a)} \cdot \BE^{(b)} e^{ 2 j \omega t }
}
\end{aligned}

Taking the difference,

\label{eqn:energyMomentumWithMagneticSources:640}
\begin{aligned}
\boldsymbol{\mathcal{E}}^{(b)} \cdot \partial_t \boldsymbol{\mathcal{E}}^{(a)}
-\boldsymbol{\mathcal{E}}^{(a)} \cdot \partial_t \boldsymbol{\mathcal{E}}^{(b)}
&=
\inv{2} \textrm{Re}
\lr{
j \omega \BE^{(a)} \cdot { \BE^{(b)} }^\conj
– j \omega \BE^{(b)} \cdot { \BE^{(a)} }^\conj
+ j \omega \BE^{(a)} \cdot \BE^{(b)} e^{ 2 j \omega t }
– j \omega \BE^{(b)} \cdot \BE^{(a)} e^{ 2 j \omega t }
} \\
&=
– \omega \textrm{Im}
\lr{
\BE^{(a)} \cdot { \BE^{(b)} }^\conj
+ \BE^{(a)} \cdot \BE^{(b)} e^{ 2 j \omega t }
},
\end{aligned}

so we have

\label{eqn:energyMomentumWithMagneticSources:660}
0
=
{
\left[
\BE^{(a)} \cross {\BH^{(b)}}^\conj
-\BE^{(b)} \cross {\BH^{(a)}}^\conj
}
+
\omega \textrm{Im}
\lr{
\epsilon_0
\BE^{(a)} \cdot { \BE^{(b)} }^\conj
+
\mu_0
\BH^{(a)} \cdot { \BH^{(b)} }^\conj
}
+ \textrm{Re}
\lr{
-\BH^{(a)} \cdot { \BM^{(b)} }^\conj
+\BH^{(b)} \cdot { \BM^{(a)} }^\conj
-\BE^{(b)} \cdot { \BJ^{(a)} }^\conj
+\BE^{(a)} \cdot { \BJ^{(b)} }^\conj
}
\right]
}_{\textrm{av}}.

Observe that the perfect cancellation of the time derivative terms only occurs when the cross product differences were those of the phasors. When those cross differences are those of the actual fields, like those in the Poynting theorem, there is a frequency dependent term is that expansion.