Most general non-relativistic Hamiltonian for hydrogenic atom in quantised electromagnetic field In many textbooks, the derivation of the energy levels in a hydrogen atom starts from the basic Hamiltonian $H = \frac{\mathbf{p}^2}{2m} + \frac{e^2}{4\pi\epsilon_0\mathbf{r}}$ and then adds relativistic fine structure or hyperfine structure as corrections to these basic energy levels. The effect of an electromagnetic field is then often included by modifying the Hamiltonian to $H = \frac{\left(\mathbf{p} - e \mathbf{A}\right)^2}{2m} + e \phi$ for electromagnetic potentials $\mathbf{A}$ and $\phi$.  
What I find strange about this approach is that from the very beginning, the interaction between proton and electron is electrostatic (thus requiring the presence of an electromagnetic field), so why is the quantised electromagnetic field only added later and not right at the start? Also, the energy of the spins of the particles in an external magnetic field does not seem to be included in this treatment.
What I wonder is the following: If we do not start from the basic Hamiltonian and add all terms later, but set out to derive from first principles the energy levels for a proton and a (bound) electron in a quantised electromagnetic field (generated by proton and electron, but also due to external radiation that might be present), what is the Hamiltonian that describes this system? The proton and electron should be considered non-relativistic, so should be regarded as quantised in the sense that observables should be operators, but without the need for quantum field theory to describe proton and electron. How could relativistic effects then be included perturbatively in this general treatment?
EDIT: Just to clarify, the situation that I am interested in is the following: Two charged, spin-carrying non-relativistic particles (proton and electron) move under the influence of an internal electric field arising from their charges and an external field due to radiation. Both particles should be included in the Hamiltonian, so the proton will also have a coupling to any electromagnetic fields. Questions: How is the internal field treated as a quantised electromagnetic field, leading to a potential energy? Which terms in the Hamiltonian does the interaction with the external field lead to? What is therefore the most general Hamiltonian for a non-relativistic proton-electron system moving in a quantised electromagnetic field?
Any references to textbooks or journal articles with a detailed treatment of the most general and (within the limits of neglecting relativistic effects) most exact Hamiltonian describing such a system would be very useful.
 A: What you want, I think, is starting from Dirac equation, find the non-relativistic limit, which is basically Pauli equation, plus spin-orbit coupling, plus the Darwin term. The classic way to do that is through the Foldy-Wouthuysen transformation [FW50]. The paper takes a bit of extra effort to read because they use old-fashioned notations for the Dirac machinery but manageable from memory. I did quote the Wikipedia pages mostly for the bibliography, as I don't find them particularly easy to learn from.
In any case, the final answer is a series in $v/c$, the classic terms of which being given by the following Hamiltonian, taken from the classic textbook by Cohen-Tannoudji, Diu and Laloë [CTDL77, chapter 12, appendix B, eqn (B-I)]:
$$H = m_ec^2 + \underbrace{\frac{P^2}{2m_e} + V}_\text{classic} \underbrace{- \frac{P^4}{8m_e^3c^2}}_\text{relativistic momentum} + \underbrace{\frac{1}{2m_e^2c^2}\frac{1}{R}\frac{\partial V}{\partial r}L\cdot S}_\text{spin orbit}+\underbrace{\frac{\hbar^2}{8m_e^2c^2}\Delta V}_\text{Darwin}$$
In the presence of an electromagnetic field $A$, the prescription is $P\to P-eA$, and then the extra term
$$-2\mu_B \frac{1}{\hbar}S\cdot B.$$
This correction has to be used to take into account the magnetic field from the magnetic moment of the proton, actually. 
Finally, there is the contribution of the electromagnetic field itself. This is a very simple term actually, just the equivalent of a sum of harmonic oscillator, one per photon mode,
$$H_{em} = \sum_i\hbar\omega_i\left(a_i^+a_i\frac{1}{2}\right).$$
So here, $a_i^+$ is a creation operator and $a_i$ is a destruction operator, both for a photon of wave vector $k_i$, energy $\omega_i$, and helicity $\epsilon_i$, perpendicular to $k_i$: this is the point of using the Colomb's gauge, it completely decouples the time and longitudinal degrees of freedom of the photons, which can then be ignored.
[FW50] Leslie L. Foldy and Siegfried A. Wouthuysen. On the Dirac theory of spin 1/2 particles and its non-relativistic limit. Phys. Rev., 78:29–36, Apr 1950.
[CTDL77] Claude Cohen-Tannoudji, Bernard Diu, and Franck Laloë. Quantum mechanics. Wiley, New York, NY, 1977. Trans. of : M ́ecanique quantique. Paris : Hermann, 1973.
A: You'd need to solve Dirac's equation
$$
(i\hbar\gamma^\mu\partial_\mu - mc)\phi = 0
$$
Here's a pretty good reference on how to do it for the Hydrogen atom
