Reference Request: Full mathematical treatment of Schrödinger evolution of hydrogen atom

Though I have not found a full reference and am still looking for one, I would like to post a partial answer to the above question:

1. The answer is a clear no. The counterexample is given in Sec. 18.3 in Hall, B. C., "Quantum Theory for Mathematicians", Vol. 267 (Springer, 2013). Basically, first observe that the $\Psi_{n l m}$ all have negative energy. Thus, any function in $L^2(\mathbb{R}^3,\mathbb{C})$ that can be expanded as a series in terms of those functions must have negative energy too. Then consider, for instance, the function $$x \mapsto \exp{\left( - \frac{(x-x_0)^2}{2} + i k \cdot x \right) }$$ with $x_0, k \in \mathbb{R}^3$. It is in $L^2(\mathbb{R}^3,\mathbb{C})$ and even in $\mathcal{S}(\mathbb{R}^3,\mathbb{C})$. Roughly speaking, this represents a wave function centered at $x_0$ traveling in the direction of $k$ at time $0$, and the magnitude of $k$ is proportional to the respective (relative) velocity. Thus, if we choose $|k|$ large enough, we can make the respective expectation value of the energy positive. Thus the function cannot be written in the above manner. Again roughly speaking, the wave function will have non-zero positive energy parts and those will scatter. As Hall writes, the $\Psi_{n l m}$ only span (in the series sense) the negative energy subspace of $L^2(\mathbb{R}^3,\mathbb{C})$. That subspace is closed by definition, and thus a Hilbert space in its own right.
2. Since the lack of square-integrability is actually due to the radial part only, we can still use the basis of spherical harmonics for the angular part. In Takhtajan's book mentioned above, the positive energy "eigenfunctions" of the radial part (up to a constant factor) are $f_{k l} (r) /r$ with energy $E(k) = \hbar^2 k^2/2 m$ and $k>0$. (The latter $m$ denotes the mass here, not the quantum number $m$; blame theoretical physicists for bad notation, not me.) The $f_{k l} (r)$ are not square integrable over $(0,\infty)$ with respect to $r^2 dr$. Asymptotically, their behavior is sinusoidal. Assuming the Gelfand triple above, and since generalized eigenfunctions should be viewed as (tempered) distributions, we $^*$should$^*$ (i.e. do not ask me for a proof) then have that any wave function $\Psi$ in the positive energy subspace of $L^2(\mathbb{R}^3,\mathbb{C})$ can be written as $$\Psi (r, \theta, \phi) = \sum_{l,m}\left( \int_{0}^\infty \operatorname{d}\negthinspace {k} \, \frac{f_{k,l}(r)}{r} \, \Phi_l(k) \right) Y_l^m(\theta, \phi) \, ,$$ where $\Phi_l \in \mathcal{S}(\mathbb{R}^3,\mathbb{C})$ for all $l$. The time evolution of $\Phi$ $^*$should$^*$ then be given by $$\Psi (t, r, \theta, \phi) = \sum_{l,m}\left( \int_{0}^\infty \operatorname{d}\negthinspace {k} \, \frac{f_{k,l}(r)}{r} \, \Phi_l(k) \, e^{- \mathfrak{i} E(k) t / \hbar}\right) Y_l^m(\theta, \phi) \, .$$ The dependence of $\Phi$ on $l$ $^*$should$^*$ come from there being an inverse relationship between any purely radial function $\tilde \Psi \in L^2(\mathbb{R}^3,\mathbb{C})$ and $\tilde{\Phi}_l \in \mathcal{S}(\mathbb{R}^3,\mathbb{C})$ whenever $$\tilde \Psi (r) = \int_{0}^\infty \operatorname{d}\negthinspace {k} \, \frac{f_{k,l}(r)}{r} \, \tilde{\Phi}_l(k) \, .$$ That is, $\tilde{\Phi}$ should depend on $l$ because $f$ does (and so should a more general $\Phi$ above).