Reference Request: Full mathematical treatment of Schrödinger evolution of hydrogen atom This is to request a reference on the above problem. I have checked various references in the mathematical physics literature and could not really find a clear, satisfactory answer to some questions, or at least I was unable to extract them myself.
In particular, I would like to know the following:

*

*It is known among mathematical physicists that most problems in quantum mechanics cannot be fully dealt with by just looking at a particular (separable) Hilbert space and are more appropriately dealt with in terms of a Gelfand triple (one reason being that $L^2$ functions are not differentiable). I suspect that in the above problem a "natural choice" for that triple is
$$\mathcal{S}(\mathbb{R}^3,\mathbb{C}) \subset L^2(\mathbb{R}^3,\mathbb{C}) \hookrightarrow \left( \mathcal{S}(\mathbb{R}^3,\mathbb{C}) \right)^* \, ,$$
where on the left we have the space of complex-valued Schwartz functions in $\mathbb{R}^3$, in the center we have the respective space of square-integrable functions, and on the right the respective tempered distributions. I base that suspicion on the fact that the common stationary states $\Psi_{n l m}$ are elements of $\mathcal{S}(\mathbb{R}^3,\mathbb{C})$. Is what is usually assumed in most physics textbooks mathematically correct, namely that any $\Phi \in L^2(\mathbb{R}^3,\mathbb{C}) $ admits a series representation of the form
$$\Phi = \sum_{n, l, m} a_{n l m} \, \Psi_{n l m} \, ?$$
If so, then the time evolution is well-defined in that problem for any such $\Phi$.

*This concerns the so-called scattering states: If one solves the time-independent Schrödinger equation for the above problem one may ask for stationary states of positive energy. Formal solutions can be found, for instance, in Takhtajan's book "Quantum Mechanics for Mathematicians". After reflecting on this question a bit, I guess the fact that those are not square-integrable means that there are no such states --- which makes physical sense, since the only way one should get a stationary state here is if the particle is somehow "trapped" by the nucleus. However, for the free particle we can solve that equation and treat the resulting states as generalized eigenfunctions in the sense that we may define
$$\Psi(x) = \frac{1}{\sqrt{2 \pi \hbar}^3} \int_{R^3} \operatorname{d}^3 p \, 
e^{\mathfrak{i}  p \cdot x / \hbar} \, \Phi(p) \, ,$$
which is, of course, just a Fourier transform. Since we can obtain $\Phi$ for given $\Psi \in L^2(\mathbb{R}^3,\mathbb{C})$ via inverse Fourier transform, the time evolution of $\Psi$ is just given by
$$\Psi(t,x) = \frac{1}{\sqrt{2 \pi \hbar}^3} \int_{R^3} \operatorname{d}^3 p \, 
e^{\mathfrak{i}  \left( p \cdot x- \omega(p) t \right) / \hbar} \, \Phi(p) \, ,$$
where $\omega(p)=p^2/2m \hbar$. Does a similar procedure work for the hydrogen atom?

EDIT: The above framing of the problem is a slight bit off. As can be found it any textbook, the ground state wave function is proportional to $e^{-|x|/a_0}$ with $a_0$ denoting the Bohr radius. Due to lack of differentiability at the origin, this function is not an element of $\mathcal{S}(\mathbb{R}^3,\mathbb{C})$. The problem seems to be fixable, however, by using $\mathcal{S}(\mathbb{R}^3 \setminus \lbrace 0 \rbrace,\mathbb{C})$ instead.
 A: 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:

*

*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.

*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).

A: The radial equation that results from separation of variables for the non-relativistic Schrodinger equation for the single electron and proton Hydrogen isotope has the form
$$
    \frac{1}{r^2}\left[-\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)+l(l+1)R-\frac{\mu\epsilon^2 r}{2\pi\epsilon_0\hbar^2}R\right]\psi=\lambda\frac{2\pi}{\hbar^2}\psi .
$$
The integer $l$ is another separation constant, which is an integer. The above may be written as
$$
        H\psi= \lambda\frac{2\mu}{\hbar^2}\psi,\;\; l=0,1,2,3,\cdots. \\
%%    \frac{1}{r^2}\left[-\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)+l(l+1)R-\frac{\mu\epsilon^2 r}{2\pi\epsilon_0\hbar^2}R\right]=\lambda\frac{2\mu}{\hbar^2}R,\;\; l=0,1,2,3,\cdots.
$$
There are square integrable solutions of this equation iff $\lambda$ is one of the following values
$$
            \lambda_n=-\frac{\mu e^4}{32\pi^2\epsilon_0^2\hbar^2n^2},\;\;\; n=l+1,l+2,l+3,\cdots.
$$
These correspond to bound states of the Hydrogen isotope, and they cluster at $0$ from below. By contrast, every $\lambda \ge 0$ is in the continuous spectrum of $H$, which corresponds to the unbound states of the system. As with the Fourier transform, which is built on eigenfunctions of $\frac{1}{i}\frac{d}{dx}$, integrals of eigenfunctions with respect to the eigenvalue parameter over any finite, positive-length interval yield square-integrable functions that are nearly eigenfunctions, if the interval of integration is small (the smaller the interval in $\lambda$, the closer the result is to an eigenfunction with eigenvalue $\lambda$.) This is analogous to the situation where the following is square integrable, and is an approximate eigenfunction of $\frac{1}{i}\frac{d}{dx}$ with approximate eigenvalue $a$:
$$
                \int_{a-\epsilon}^{a+\epsilon}e^{i\lambda x}d\lambda
$$
The discrete, bound states, together with the continuum of unbound states combine to give a complete $L^2$ expansion for functions in $L^2(\mathbb{R})$. These expansions require both discrete sets of eigenfunctions corresponding to bound states, as well as the non-square-integrable classical eigenfunctions in order to be complete.
