Why maximum energy transfer at natural frequency even if max amplitude occurs below $f_0$
Interesting question. Note that the power transfer is $$P=\frac{1}{\tau}\int_0^\tau F(t)v(t)\,\mathrm dt$$ where $F(t)$ is the applied force, $v(t)$ is the velocity of the oscillation, and $\tau$ is the length of the period of the oscillation (this is just a simple extension of the usual rule for work, $w=Fd$, to time-varying systems).
The usual forced oscillator $$m x''(t)=-kx(t)-\gamma x'(t)+A\sin(\omega t)$$ has steady-state solution $$x(t)=\frac{A }{\sqrt{\left(k-m\omega^2\right)^2+\gamma ^2 \omega ^2}}\sin \left(\omega t+\tan ^{-1}\left(k-m \omega ^2,-\gamma \omega \right) \right).$$
Since the user-applied force is $F(t)=A\sin(\omega t)$, the velocity is $v(t)=x'(t)$ and $\tau=2\pi/\omega$, we obtain \begin{align}P&=\frac{\omega}{2\pi}\int_0^{2\pi/\omega} A\sin(\omega t)\frac{A \omega \cos \left(\omega t+\tan ^{-1}\left(k-m \omega ^2,-\gamma \omega \right)\right)}{\sqrt{\left(k-m\omega^2\right)^2+\gamma ^2 \omega ^2}}\,\mathrm dt\\&=\frac{\frac{1}{2}A^2 \gamma \omega ^2}{\left(k-m \omega ^2\right)^2+\gamma ^2 \omega ^2}\,.\end{align} Solving $$\frac{\mathrm dP}{\mathrm d\omega}=0$$ for $\omega$ yields $$\omega=\sqrt{\frac{k}{m}}$$ which is exactly the natural frequency.
Thus, somewhat paradoxically, maximum power transfer occurs when the forcing is at the natural frequency, even though the vibrational amplitude is not maximum.
(Minor note: I am using the two-argument arctangent function $\tan^{-1}(a,b)$ as it is defined in Mathematica, ArcTan[a,b]
).