How (if) are the hydrogen's "external" and "internal" wavefunction connected? We can assign to a hydrogen atom an "internal" wavefunction:
$$\psi(r,\theta,\varphi)=R(r)Y(\theta,\varphi)$$
We also assign an "external" wavefunction, describing the behavior of the atom moving freely through space. The atom is viewed as a point particle in this case though:
$$ \psi(\mathbf{r}, t) = Ae^{i(\mathbf{k}\cdot\mathbf{r}-\omega t)} $$
The wavefunction of a localized particle (in this case, hydrogen) is made out of a distribution of different wavelength plane wave wavefunctions. A Gaussian wave packet.
Are these two somehow connected when constructing the total wavefunction of the atom? Can the Gaussian wavepacket (the external wavefunction) and "internal" wavefunction be somehow connected to obtain some overall total wavefunction?. Or is this distinction superficial, i.e., can both be used only in separated situations?
 A: A wavefunction is the mathematical solution of wave equations, and there are more than one equations used in building quantum mechanical models for atoms etc.
What you call "internal wave function" for hydrogen , is the solution of the   specific equation with the Coulomb $1/r$ potential for the hydrogen atom, i.e. for the masses and charges  of electron and proton .
What you call an "external" wavefunction is a solution of the free wave equation: no potential so the simplest solution is the plane wave. This might  be used as a representative of a quantum size particle except, there is equal probability to be anywhere in space time,  so it cannot be used  to represent a real particle at (x,y,z,t). One has to use more complicated solutions , the wave packet solution, in order to really represent a localized particle, i.e.  if one needs to represent a neutral hydrogen atom.
(In quantum field theory, plane waves of each particle in the table  ( electron, quark  ...)constitute the field on which creation and annihilation operators work , in the mathematics of the Feynman diagram that allow to compute particle interactions.)
A: In principle you can solve the wave equation in the atomic rest frame, using the Born-Oppenheimer approximation to keep things tractable. Then you can apply a Galilei or Lorentz transformation to get a moving atom. However the Schrödinger equation is not covariant under a Galilei transformation. You perhaps can use the Schrödinger group. However, see my arxiv paper which criticizes this approach. For the Dirac equation the approach should give you you the wave function that you are looking for. I doubt that it can be separated in the way that you propose.
A: Short answer:  yes, they're equally real, and they combine to form the full wavefunction for the system.
The full wavefunction $\Psi(\mathbf{r}_e, \mathbf{r}_p)$ for an electron and proton interacting with each other is the solution to the 6-variable PDE
$$
- \frac{\hbar^2}{2 m_e} \nabla^2_{r_e} \Psi - \frac{\hbar^2}{2 m_p} \nabla^2_{r_e} \Psi - \frac{k e^2}{|\mathbf{r}_e - \mathbf{r}_p|}\Psi = E \Psi,
$$
where $m_e$ and $m_p$ are the respective masses of the electron and proton, $\mathbf{r}_e$ and $\mathbf{r}_p$ are their respective positions, and $\nabla^2_{r_e}$ stands for the Laplacian involving derivatives with respect to $\mathbf{r}_e$ (and similarly for $\nabla^2_{r_p}$.)
Trying to solve this equation as it stands is pretty much impossible.  There are six different variables in play (three for each particle) and all of the equations are coupled together via the potential term.  However, a clever change of coordinates can make our lives simpler.  We can define
\begin{align*}
\mathbf{R} &= \frac{m_e \mathbf{r}_e + m_p \mathbf{r}_p}{M} = \text{position of the center of mass} \\
\mathbf{r} &= \mathbf{r}_e - \mathbf{r}_p = \text{separation vector between electron & proton}
\end{align*}
where $M = m_e + m_p$.  If we rewrite the Schrödinger equation in terms of these variables, it becomes
$$
- \frac{\hbar^2}{2 M} \nabla^2_{R} \Psi - \frac{\hbar^2}{2 \mu} \nabla^2_{r} \Psi - \frac{k e^2}{r}\Psi = E \Psi,
$$
where $\mu$ is the reduced mass defined via $\mu^{-1} = m_e^{-1} + m_p^{-1}$.
Why does this help us?  The PDE is separable now!  Specifically, we can look for a solution of the form
$$
\Psi(\mathbf{r}_e, \mathbf{r}_p) = \psi(\mathbf{R}) \phi(\mathbf{r})
$$
where the functions $\psi$ and $\phi$ satisfy
$$
- \frac{\hbar^2}{2 M} \nabla^2_{R} \psi = E_R \psi \\
- \frac{\hbar^2}{2 \mu} \nabla^2_{r} \phi - \frac{k e^2}{r}\phi = E_r \phi
$$
If we can find these solutions, then their product will satisfy the full 6-D Schrödinger equation with $E = E_R + E_r$.  But the solutions to the first equation are just plane waves in $\mathbf{R}$, while the solutions to the second equation are the familiar solutions for the hydrogen atom.  So the energy eigenstates of the system would be of the form
$$
\Psi(\mathbf{r}_e, \mathbf{r}_p) = e^{i \mathbf{k} \cdot \mathbf{R}} R(r) Y(\theta,\phi)
$$
with $\mathbf{R}$ and $\mathbf{r}$ defined in terms of $\mathbf{r}_e$ and $\mathbf{r}_p$ as above.
You could then construct a Gaussian wavepacket by taking a superposition of several such plane waves in $\mathbf{R}$.  The resulting wavefunction could then be used to predict (say) the expected position of the electron, by calculating the integral
$$
\langle \mathbf{r}_e \rangle = \int \mathbf{r}_e \Psi^*(\mathbf{r}_e, \mathbf{r}_p) \Psi(\mathbf{r}_e, \mathbf{r}_p) \, d^3 \mathbf{r}_e \,d^3 \mathbf{r}_p
$$
Note that both the "plane-wave" and "hydrogenic" parts of $\Psi$ would need to be taken into account, since they both depend on $\mathbf{r}_e$ via their dependence on $\mathbf{R}$ or $\mathbf{r}$.
