Hamiltonian for charge particle in electromagnetic field in dipole approximation Why can the Hamiltonian
$$
H = \frac{1}{2m}\left(\vec{p}-q\vec{A}\right)^2+q\Phi + V
$$
be transformed to
$$
H = \frac{1}{2m}\vec{p}^2 - q \vec r \vec E + V
$$
in the dipole approximation, in which the vector potential $\vec{A}$ (and, as a consequence, the electric field $\vec{E}$) are assumed to be constant in space?
So far, I have found three approaches to derive the second Hamiltonian.

*

*Expand the quadratic expression and neglect $q^2 \vec{A}^2$. While this is applied frequently in literature, it is actually not required in dipole approximation.

*Transform the resulting wave functions $\psi' = \exp(-\mathrm{i}\hbar^{-1}q\vec{r}\vec{A})\psi$, as described, e.g., in Meystre "Elements of Quantum Optics". But why is such a transform valid? Clearly, the probability distribution $|\psi'|^2 = |\psi|^2$ does not change, but it is not exactly a gauge transform, is it?

*Maria Goeppert-Mayer [1] states that the Lagrangian corresponding to the first Hamiltonian in this post is equivalent to another Lagrangian, which is derived by adding the total differential $\mathrm{d}_t(\vec{r} \vec{A})$. The second Lagrangian can be converted to the second Hamiltonian in this post. But why is this addition valid?

[1] Göppert‐Mayer, Maria. "Über Elementarakte mit zwei Quantensprüngen." Annalen der Physik 401.3 (1931): 273-294.
 A: In the hope that it may be helpful to some, I would like to post a short summary of my findings.
We note that the canonical momentum $\vec{p}'$ is different from the mechanical momentum $\vec{p} = m \vec{v}$ in the first Hamiltonian. We aim to derive a form in which both momenta coincide. This can be achieved by adding a total time differential to the Lagrangian, as described in the paper by M. Goeppert-Mayer. In this post and the references therein, it is explained thoroughly why the Euler-Lagrange equation is invariant to this operation. Thanks to @Philip for pointing it out to me.
In the scope of the Hamiltonian formalism, this operation corresponds to a canonical transform. Here, I found this answer very helpful.
Basing on this answer, we use the transform
$$\vec{r}' = \vec{r}$$
$$\vec{p}' = \vec{p} + \nabla_\vec{r} F $$
$$H' = H - \partial_t F$$
and define $F = - q \vec{r} \vec{A}$.
We note that the vector potential $\vec{A}$ is assumed to be constant in space (at least in the dimensions of the quantum mechanical system), so that
$$\vec{p}' = \vec{p} - q A$$
and write for the new Hamiltonian
$$H'(\vec{r}, \vec{p}', t) = \frac{1}{2m}\vec{p}'^2 + V + \partial_t q \vec{r} \vec{A} = \frac{1}{2m}\vec{p}'^2 + V - q \vec{r} \vec{E}.$$
Note 1: The scalar potential $\varphi$ is usually assumed to be zero in the derivations I have seen so far.
Note 2: When addressing the subtlety mentioned here it is helpful to know that $\partial_{\vec{r}}^2 F = 0$.
