Good question. It exemplified a key issue I think that exists between what (particle) physicists typically do, and what mathematicians would consider well-defined. The issue comes about due to the fact that the system in question is infinite dimensional.
To illustrate the issue, let's generically say that $\hat{H}$ is a free-Hamiltonian with some associated equation of motion $Df=0$. (One can use the Klein-Gordon system as an example.) The eigenstates of the Hamiltonian $\phi$ are typically functions that cannot be normalized. However, being eigenstates of a Hermitian operator, they do obey some orthogonality condition. There is expected to be a very large degeneracy for such an infinite dimensional system. However, a suitable way to define orthogonal eigenstates within this degeneracy is to use plane wave eigenfunctions $\phi(\mathbf{k})$ which are parameterized by the wavevector $\mathbf{k}$. Now we can write the orthogonality condition as
$$ \langle\phi(\mathbf{k}_1),\phi(\mathbf{k}_2)\rangle = 2\pi\delta(\mathbf{k}_1-\mathbf{k}_2) . \tag{1} $$
Obviously these functions are not normalizable.
To resolve this issue, we can look at the solutions of the equation of motion instead of the eigenfunctions of the Hamiltonian. While all eigenfunctions are solutions of the equation of motion, the converse is not in general true. Therefore, there are many more solutions than eigenfunctions. One can select those solutions that can be normalized. The way to do that is to start with a space of well-defined normalizable functions $F(\mathbf{k})$ (for example Schwartz space), which we call spectra. While the plane waves also depend on the spacetime coordinates $\phi(\mathbf{x},t;\mathbf{k})$, these spectra only depend on the wavevector. So the space of solutions that we are looking for can now be defined by
$$ f(\mathbf{x},t) = \int F(\mathbf{k}) \phi(\mathbf{x},t;\mathbf{k}) d^3k . \tag{2} $$
This is somewhat related to an inverse Fourier transform, but the fact that we increase the number of dimensions need to be taking into account. It follows from the fact that the plane wave satisfy the equation of motion, which introduces a dimensional reduction. The space of solutions $f(\mathbf{x},t)$ inherits the properties of the space of spectra. The plane waves only serve to provide the orthogonality and completeness properties.
The space of solutions is an inner product space with an inner product inherited from the inner product between plane waves. For instance, consider an equal-time inner product between solutions:
$$ \begin{aligned}
\langle f_1, f_2 \rangle_t & = \int F_1^*(\mathbf{k})
\phi^*(\mathbf{x},t;\mathbf{k}) F_2(\mathbf{k}') \phi(\mathbf{x},t;\mathbf{k}') d^3k' d^3k d^3x \\
& = \int F_1^*(\mathbf{k}) \delta(\mathbf{k}-\mathbf{k}')
F_2(\mathbf{k}') d^3k' d^3k \\
& = \int F_1^*(\mathbf{k}) F_2(\mathbf{k}) d^3k .
\end{aligned} \tag{3} $$
The properties of the spectra (Schwartz functions) allow us to change the order of integration. The resulting inner product is now simply defined as an integral over the spectra of the two solutions with one be complex conjugated. If we consider a complete space for the spectra (for example $L^2(\mathbb{R}^3,\mathbb{C})$ with Schwartz space being a dense subset), the completeness of the space of solutions thus follows from the properties of Fourier transforms, being an isometry between Hilbert spaces.
Response to comments: From the discussion in the comments it follows that the definition of the space of spectra is not a problem. Just to summarize, we can nominally take it as $L^2(\mathbb{R}^3,\mathbb{C})$. Eventually we'll want to restrict it to $\mathcal{S}(\mathbb{R}^3)$. Note however, that at this point the idea of Lorentz invariance has no meaning. The relevant measure is just the standard Lebesgue measure.
The question now is how we define the "space of solutions." Nominally the solutions of the equation (or the eigenfunctions of the Hamiltonian) can be specified as the set of all plane waves $$\phi(\mathbf{x},t;\mathbf{k}) = \exp(-i\omega t + i \mathbf{k}\cdot\mathbf{x}) , $$
where
$$ \omega = c\sqrt{k_m^2+\mathbf{k}\cdot\mathbf{k}} , $$
with $k_m$ being the wavenumber associated with the mass, is the dispersion relation (or mass-shell condition). As we already stated above, these plane waves obey an orthogonality condition given in (1). These plane waves also allow us to introduce the notion of Lorentz invariance.
Now we can define the solutions as in (2). (Let's represent the space of all such solutions by $\mathcal{F}$.) What does it mean? One may want to consider such solutions as elements of $L^2(\mathbb{R}^4,\mathbb{C})$, but not all functions in $L^2(\mathbb{R}^4,\mathbb{C})$ are solutions that can be expressed as (2). Therefore, the Hilbert space $\mathcal{H}^4\equiv L^2(\mathbb{R}^4,\mathbb{C})$ does not represent $\mathcal{F}$. Nor can we say that $\mathcal{F} \subset \mathcal{H}^4$, because the inner product is not the same. Instead the inner product is defined in (3) with the aid of (1). Therefore, $\mathcal{F}$ is isomorphic to $\mathcal{H}^3\equiv L^2(\mathbb{R}^3,\mathbb{C})$, considered as the closure of $\mathcal{S}(\mathbb{R}^3)$, with the isomorphism defined via (3).
So what is different about $\mathcal{F}$? The solutions in $\mathcal{F}$ allow us to introduce Lorentz transformations. Note however, that the definition of the inner product in (3) is not Lorentz invariant.
Did I manage to address all the issues now?
