Born Oppenheimer Approximation: Why can any molecular state be represented as a linear combination of electronic states? in the Born Oppenheimer Approximation, one expands the molecular wavefunction $\Psi(x,X)$ in terms of the electronic wavefunctions $\phi(x;X)$:
$$\Psi(x,X)= \sum_k(c(X)_k\phi(x;X)_k)$$
($x$ are the electronic coordinates and $X$ are the nucleonic coordinates)
Now, since the electronic wavefunctions are eigenstates of the electronic Hamiltonian, the constitute a complete basis of the electronic space. Thus any electronic wavefunction can be expanded in terms of the eigenfunctions.
But, how can we be sure that any molecular wavefunction can be expanded in terms of the electronic wavefunctions? How can we be sure that the molecular Hilbert space is not larger than the space which is spanned by the eigenstates of the electronic Hamiltonian?
 A: I know this question was asked a long time ago, but since I thought very hard about the same question today and didn't find the other answer very helpful, I decided to write my own. The problem is that OP is not really asking about why the product form used in the BO approximation is valid, but how the given expansion (which is claimed to be exact, see e.g. the scan of a book chapter provided in this post: Non Adiabatic Coupling Term in Born Oppenheimer Approximation) is justified.
Even though this problem confused me a lot, it is actually quite simple, and I understood it best using some mathematical language. I call the space of electronic coordinates $A$ (so $x \in A$) and $\Psi(x, X)$ is the exact solution of the complete (electronic and nuclear) Schrödinger equation. Then for each set of nuclear coordinates $X$, we can define a function
\begin{equation}
\Psi_X: A \rightarrow \mathbb{C}, \quad \Psi_X(x) := \Psi(x, X).
\end{equation}
We know that for all $X$, the set of eigenfunctions of the electronic Hamiltonian, {$\phi(x;X)_k$}, is a complete basis for electronic wave functions. We now expand each electronic wave function $\Psi_X$ in a different set of basis functions, namely
\begin{equation}
\Psi_X(x) = \sum_k c(X)_k \phi(x;X)_k,
\end{equation}
where $c(X)_k$ is the expansion coefficient belonging to the $k$th basis function associated with the electronic wave function parametrized by $X$. But since $\Psi_X(x) = \Psi(x, X)$ we already obtained what OP asked for.
A: The key word here is approximation.  Its not the that the hilbert space of the full molecule is not larger, its that the velocities of the much more massive nucleons can be neglected and only their positions need to be considered.  In fact, the BO approximation breaks down in some cases when non-adiabatic effects are not negligible (e.g. loss or gain of energy due to changes in electronic orbits).  An important point that Born and Oppenheimer make in their paper is their approximation is possible because earlier classically based approximations overstated the nuclear contribution because it neglected spin, which is purely a quantum phenomenon with no classical analog.  This causes the nuclear contribution to be of the fourth order instead of the second order.  This is another excellent case to show where classical mechanics breaks down in the description of the universe.
I guess to answer your question though, the hilbert space of the full system is actually larger then the electronic states, its just that under certain conditions some nuclear degrees of freedom can be ignored.
