We can solve the Schrodinger equation for the Hamiltonian operator from the classical Hamiltonian of hydrogen bound state, consisting of proton and electron attracting each other electrodynamically, to get the eigenfunction, which is corresponding to stationary bound state we wanted to find.
But It must also be possibly explained through QED, but maybe rather complicated way. I can't even imagine the way it does. Could anyone explain it for me possibly qualitatively?
Of course in the spectrum of QED there is hydrogen atom. The problem is what to compute in the case of bound states and how to compute. The usual approach to a general QFT is perturbative, meaning that you start from free fields (in this case electron and photon) and then you think the compulg coupling between these fields (in this case $\bar\psi A_\mu\gamma^\mu\psi$) as "little" (all this statement can be made mathematically rigorous using spectral theory...). Intuitively this means that you are thinking that the particle described by fields are "almost free". This is the case in a scattering process where you assume that the initial and final state are "free particle" but, of course it cannot be true for a bound state. So all the technique you probably learned from perturbation theory (Faynman graph and so on) cannot be applied to this case. You must focus of some nonperturbative aspect of the theory. If you where able to compute the spectral density (the one that appear in Källén–Lehmann) representation you would notice a delta function in correspondence of the hydrogen atom energy. Then other signals of the bound states can be seen in some pole of scattering matrix. But this are hard or impossible to compute. An approach to this problem is the Bethe-Salpeter equation.
