Solution (source code)

= Solution

For the <isotropic harmonic oscillator flow>, <Hamilton's equations> are $\dot X=V$, $\dot V=-\omega^2X$. In dimension three these are three identical uncoupled pairs. For $\omega\ne0$, put $c_t=\cos(\omega t)$, $s_t=\sin(\omega t)$. The solution from $(x,v)$ at time zero is
$$
\boxed{X(t)=c_tx+\frac{s_t}{\omega}v,\qquad
V(t)=-\omega s_tx+c_tv.}
$$
The inverse flow is obtained by replacing $t$ by $-t$:
$$
S_{-t}(x,v)=\left(c_tx-\frac{s_t}{\omega}v,\ \omega s_tx+c_tv\right).
$$
Integrating the source along the backward characteristic gives the <Duhamel formula for Hamiltonian transport>:
$$
\boxed{f(t,x,v)=
f_0\!\left(c_tx-\frac{s_t}{\omega}v,\ \omega s_tx+c_tv\right)
+\int_0^t h\!\left(s,\ c_{t-s}x-\frac{s_{t-s}}{\omega}v,\
\omega s_{t-s}x+c_{t-s}v\right)\,ds.}
$$
Indeed $f(t,S_tz)=f_0(z)+\int_0^th(s,S_sz)\,ds$, and setting $z=S_{-t}(x,v)$ gives the formula. It has the prescribed initial value and differentiation along the characteristic gives the source.

The zero-frequency limit has $s_t/\omega\to t$, $\omega s_t\to0$, so $S_t(x,v)=(x+tv,v)$ and
$$
\boxed{f(t,x,v)=f_0(x-tv,v)+\int_0^th(s,x-v(t-s),v)\,ds\qquad(\omega=0).}
$$