Solution (source code)

= Solution

Define the unperturbed orbital <angular frequency> vector by
$$
\boxed{\boldsymbol\Omega(\mathbf J)=\frac{\partial H_0}{\partial\mathbf J}.}
$$
It is the rate of advance of the <angle variables> along the unperturbed <Hamiltonian flow>; the <action variables> are constant. Expand the <Collisionless Boltzmann equation> to first order in the perturbations. Since $f_0$ and $H_0$ depend only on the <action variables>, $[f_0,H_0]=0$, and
$$
\partial_t\delta f+[\delta f,H_0]+[f_0,\delta\Phi]=0,
\qquad
\partial_t\delta f+\boldsymbol\Omega\cdot\partial_{\boldsymbol\theta}\delta f
-\partial_{\mathbf J}f_0\cdot\partial_{\boldsymbol\theta}\delta\Phi=0.
$$
The omitted <Poisson bracket> $[\delta f,\delta\Phi]$ is quadratic in the perturbations. Use the <Fourier series> convention
$$
\delta f=\sum_{\mathbf k\in\mathbb Z^3}\delta f_{\mathbf k}(\mathbf J,t)e^{i\mathbf k\cdot\boldsymbol\theta},\qquad
\delta\Phi=\sum_{\mathbf k\in\mathbb Z^3}\delta\Phi_{\mathbf k}(\mathbf J,t)e^{i\mathbf k\cdot\boldsymbol\theta}.
$$
For each harmonic the linear <ordinary differential equation> is
$$
\partial_t\delta f_{\mathbf k}+i\mathbf k\cdot\boldsymbol\Omega\,\delta f_{\mathbf k}
=i\mathbf k\cdot\partial_{\mathbf J}f_0\,\delta\Phi_{\mathbf k}.
$$
Multiplication by the <integrating factor> $e^{i\mathbf k\cdot\boldsymbol\Omega t}$ and use of $\delta f_{\mathbf k}(0)=0$ give the causal <linear response in action-angle variables>:
$$
\boxed{\delta f_{\mathbf k}(\mathbf J,t)
=i\mathbf k\cdot\partial_{\mathbf J}f_0\int_0^t\delta\Phi_{\mathbf k}(\mathbf J,t')e^{-i\mathbf k\cdot\boldsymbol\Omega(t-t')}\,dt'.}
$$
The zero harmonic stays zero at this order. This formula retains the entire forcing history; replacing it immediately by a resonant <Dirac delta function> would discard the transient response.