Justifying separation of variables in solving the Schrödinger equation in 3D

Question

Consider a 3D infinite square "box", satisfying the time-independent S.E. $$ \begin{align*} -\frac{\hbar^2}{2m}\nabla^2 \psi &= E\psi, \;\;\; x,y,z\in[0,a], \\ \psi &= 0, \;\;\;\;\;\; \text{otherwise}. \end{align*} $$ The standard method for finding the allowed energies seems to be separation of variables. We write $\psi(\mathbf r) = X(x)Y(y)Z(z)$, so $$ -\frac{\hbar^2}{2m}\left(YZ\frac{\mathrm d^2 X}{\mathrm d x^2} + XZ\frac{\mathrm d^2 Y}{\mathrm d y^2} + XY\frac{\mathrm d^2 Z}{\mathrm d z^2} \right) = EXYZ, $$ and dividing by $XYZ$ on both sides yields $$ -\frac{\hbar^2}{2m}\left(\frac 1X\frac{\mathrm d^2 X}{\mathrm d x^2} + \frac 1Y\frac{\mathrm d^2 Y}{\mathrm d y^2} + \frac 1Z\frac{\mathrm d^2 Z}{\mathrm d z^2} \right) = E. $$ At this point the problem is reduced to three 1D infinite square well equations, and we get the allowed energies $$ E = \frac{\pi^2 \hbar^2}{2ma^2}(n_x^2 + n_y^2 + n_z^2), \;\;\; n_x,n_y,n_z \in \mathbb N, $$ and eigenfunctions $$ \psi = X_{n_x} Y_{n_y} Z_{n_z} = \left(\frac 2a\right)^{3/2}\sin\left(\frac{n\pi y}{a} \right)\sin\left(\frac{n\pi x}{a} \right)\sin\left(\frac{n\pi z}{a} \right). $$ My question is, how can we be sure that we've obtained all of the possible eigenvalues/eigenfunctions, and that there isn't some rogue eigenfunction/eigenvalue roaming somewhere out there, that cannot be expressed as a separable solution $\psi = XYZ$?

For the 1D infinite square well, for example, we were able to obtain the general solution to the differential equation $$ -\frac{\hbar^2}{2m}\frac{\mathrm d^2 \psi}{\mathrm d x^2} = E\psi $$ and apply boundary conditions to get the allowed energies. Hence, I can see why we must have gotten all of the possible eigenvalues and eigenfunctions in that case. But for the 3D case $\psi = XYZ$ does not seem to be the general solution, though I could be mistaken.

More importantly, separation of variables is used in solving the hydrogen S.E., by assuming that $\psi(\mathbf r) = R(r)Y(\theta,\phi)$. I have a hard time convincing myself mathematically that there is not some non-separable eigenfunction out there, with some energy that cannot be written as $-13.6 n^{-2} \; \mathrm{eV}$, though experimental results clearly suggest otherwise.

Edit: To clarify my question, I understand how $\hat H \psi = E\psi$ is derived using separation of variables from its time-dependent counterpart $\hat H \Psi = i\hbar \partial _t \Psi$, and how the solutions $\{\psi\}$ of $\hat H \psi = E\psi$ are a complete orthonormal set that can be linearly combined to construct any $\Psi$ in Hilbert space that satisfies $\hat H \Psi = i\hbar \partial _t \Psi$. ($\hat H$ is an observable operator, and QM assumes that the eigenfunctions of all observable operators are a complete set.)

My question is, in applying separation of variables to $\hat H \psi = E\psi$ again, how can we be sure that we didn't miss out on any of those $\psi$'s which make up that complete orthonormal set with which we construct solutions to $\hat H \Psi = i\hbar \partial _t \Psi$, as we do not obtain the general solution to $\hat H \psi = E\psi$, as we usually do in 1D cases?

Answer 1

In short, the set of eigenfunctions you've found is a complete orthonormal basis for the full Hilbert space of the system. That means that any 'rogue eigenfunction/eigenvalue roaming somewhere out there' must be a linear combination of the eigenfunctions you've already found.

Answer 2

It can be proven that any reasonable function can be expanded in terms of the $X$, $Y$, $Z$'s, i.e., they are complete. One can with relative ease convince oneself of this fact starting from the fact that $X_n(x)$ is complete in the function space of 1D functions, i.e., any single-variable function $f(x)$ can be expanded in terms of $X_n(x)$, which one knows from the mathematical knowledge of Fourier series.

Due to the above completeness with respect to 3D functions, there cannot be any other eigenfunctions. Because if there is another eigenfunction, this eigenfunction, which is also a 3D function, must be orthogonal to all the eigenfunctions already found and thus cannot be expanded into those eigenfunctions $X$, $Y$ and $Z$'s. This is contradictory to the above completeness which says any 3D function can be expanded in the eigenfunctions $X$, $Y$ and $Z$. Thus one knows all eigenfunctions have been found.

Answer 3

To be a little more explicit than @Emilio Pisanty.

When you write $n^2_i$ instead of $-\frac{\partial^2 }{\partial x_i^2}$, you have implicitly written the wave function as its Fourier series. The set of all the functions of this Fourier series is a complete basis (set of base functions), i.e. you can approximate every generic function as a combination of these basis functions.

If your wave function is written as a linear composition of many base functions (eigenvectors) you can't associate a specific energy value to your function, but instead a probability associated to the contribution of its projection on the base functions.