The Darwin term can be obtained from the low energy approximation ($|p|^2/m^2<<1$) of the Dirac equation of the electron in a central field. An elegant way to perform this task is by means of the Foldy-Wouthuysen Transformation. The same approximation leads also to the other terms detected in the Hydrogen atom fine structure including the spin-orbit term.
However, the relation to the zitterbewegung is also correct, at least in the heuristic level. Please see a derivation in the following work by Klaus Capelle
(Equation 4). The explanation is that due to the zitterbewegung high frequency motion, the central potential gets smeared on the average and acquires the extra Darwin term.
Charles Darwin together Gordon, were the first to compute the exact energy levels of the Hydrogen atom using the Dirac equation in 1928, just two years after Dirac discovered his equation. Darwin also performed the low energy approximation to understand the differences from the nonrelativistic treatment.
A little piece of history about the Darwin term can be found in footnote 7 of the following lecture notes which is quoted in the following:
"Sir Charles Galton Darwin (1987-1962), a British physicist, was a grandson of the Charles Darwin of evolution fame. Sir Charles was the first, with Gordon, to work out the exact energy levels of Hydrogen according to the Dirac equation and thereby discovered the eponymous term in the levels. He worked out the Lagrangian and Hamiltonian for classical motion of several interacting charges correct to O($v^2/c^2$). He also worked on statistical mechanics (Darwin-Fowler method). Later in life he took part in the Manhattan project."
Darwin's work was referred to in the original Foldy-Wouthuysen article: Phys Rev. 78 no.1, 29-36, 1950.
Best Answer
You can find the complete derivation in ref.1, by using the Dirac equation. You may want to complement it by having a look at ref.2. first, where the Dirac equation is derived from the principles of quantum electrodynamics (which is a more fundamental theory), thus obtaining a more complete picture.
Sketch: ref.2. takes the QED Lagrangian, with operator $\hat\psi(x)$, and derives the Dirac equation for the field $\psi_n(x):=\langle 0|\hat\psi(x)|n\rangle$, where $|n\rangle$ is a complete set of state-vectors and $|0\rangle$ is a vacuum state. Given the Dirac equation for $\psi_n(x)$, ref.1. performs the so-called Foldy-Wouthuysen transformation, and then takes the non-relativistic limit, thus obtaining the first-order corrected Schrödinger equation. All in all, refs.1,2 provide a first-principles derivation of the Hamiltonian OP is after (up to factors of $Z$, which are easy to reintroduce).
References.
Itzykson and Zuber, QFT, §2-2-4.
Weinberg, QFT, Vol.I, §14.1.