Eigenvalue equation and separable solutions Recently, studying Quantum Mechanics I found a doubt regarding separable solutions and eigenvalue equations for differential operators. Suppose we are considering some space $\mathcal{H}$ of functions in $\mathbb{R}^3$ with values in $\mathbb{C}$, like $L^2(\mathbb{R}^3)$. Now consider a certain differential operator $A$ in $\mathcal{H}$ and the eigenvalue equation for $A$:
$$A\psi = \lambda \psi.$$
This is a partial differential equation, if $A$ is a differential operator. Now, suppose for simplicity we already know the spectrum of $A$ and we want to solve it to find the eigenvectors.
Since it is a PDE the first strategy would be separation of variables. In that case we consider, using for example cartesian coordinates that $\psi(x,y,z)=X(x)Y(y)Z(z)$. If everything works, what we will achieve is that we will find that if $\psi$ is separable, then it is a solution with $X$,$Y$ and $Z$ given by something we found.
This, however, doesn't say the converse: that if $\psi$ is an eigenfunction of $A$, it is necessarily of the form $\psi(x,y,z)=X(x)Y(y)Z(z)$. Indeed, if $\psi_1$ and $\psi_2$ are separable solutions their sum, which is not a separable solution in general, also satisfies $A(\psi_1+\psi_2)=\lambda(\psi_1+\psi_2)$.
My question is: if $\psi$ is an eigenfunction of a differential operator $A$ for a known eigenvalue $\lambda$, is $\psi$ necessarily separable?
My whole point is: when we pick the eigenvalue equation and use separation of variables we show that "being separable implies being an eigenfunction", but the converse doesn't seem to be imediate. Is the converse true?
 A: Separation of variables works on regions that are rectangular in some particular coordinate system. The underlying operator must be separable in that coordinate system. For the Laplacian, that generally means an orthogonal coordinate system such as spherical coordinates, cylindrical coordinates, elliptic coordinates, etc.. A rectangular region in spherical coordinates can be a sphere, a spherical thick shell, a wedge, etc.. There are a couple of dozen such coordinate systems. You can add various potentials, provided their dependence is also separable (usually meaning that the it depends on only one of the coordinate variables.)
So the number of configurations where you can separate variables is limited, but it includes an important class of problems. Once you are able to separate variables, the result is an ODE in each coordinate on an interval in the corresponding variable, which is why the region where you're solving needs to be a rectangular box in the chosen coordinate system. Sturm and Liouville did a nice job in the early 1800's of characterizing the ODEs coming out of separation of variables eigenvalue problems. These are the Sturm-Liouville eigenfunction equations with eigenvalue $\lambda$ on an interval $[a,b]$:
$$
    \frac{1}{w}\left[\frac{d}{dx}\left(p(x)\frac{df}{dx}\right)+q(x) \right]=\lambda f, \\
         \cos\alpha f(a)+\sin\alpha f'(a) = 0,\\
         \cos\beta f(b)+\sin\beta f'(b) = 0.
$$
The function $w$ determines a weight for the space $L^2_w(a,b)$ where the problem is properly posed as a selfadjoint one. The inner product on this space involves this weight function, and as given by
$$
         (f,g)_w=\int_{a}^{b}f(t)\overline{g(t)}w(t)dt.
$$
The function $p$ generally comes out of a scale factor associated with the coordinate change and is positive on $(a,b)$; and $q$ is a potential. Infinite and semi-infinite regions may occur (i.e., where $a=-\infty$ and/or $b=\infty$.) If $p$ vanishes at an endpoint of $(a,b)$, then the problem is singular, and there may or may not be some type of endpoint condition at the singular endpoint.
If you're working on a finite region in all coordinates, and $p$, $w$ are strictly positive on $[a,b]$, then you have a regular Sturm-Liouville problem. In that case there exists a sequence of eigenvalues that are bounded below, and a set of mutually orthogonal eigenfunctions. The normalized eigenfunctions form a complete orthonormal basis of $L^2_{w}(a,b)$, which means that every $f \in L^2_w(a,b)$ can be expanded in a series of such eigenfunctions. An equivalent of having a complete orthonormal basis $\{ \varphi_{n} \}_{n=1}^{\infty}$ is that, if $f \in L^2_w(a,b)$, then
$$
            \int_{a}^{b}f(t)\varphi_{n}(t)w(t)dt = 0,\;\; 1 \le n \le \infty \implies f = 0.
$$
Suppose you have a 2-d region. Then you are working on $[a,b]\times [c,d]$, and you have a complete basis in each direction, say $\{ \varphi_{n}(u) \}_{n=1}^{\infty}$ and $\{ \psi_{n}(v)\}_{n=1}^{\infty}$. Then $\{ \varphi_{n}(u)\psi_{m}(v) \}_{n,m=1}^{\infty}$ is general a complete basis because
$$
       \int_{c}^{d}\int_{a}^{b}f(u,v)\varphi_{n}(u)\psi_{m}(v)w_1(u)w_2(v)dudv=0,\;\;\; n=1,2,3,\cdots,\; m=1,2,3,\cdots \\
  \implies \int_{c}^{d}\psi_{m}(v)\left[\int_{a}^{b}f(u,v)\varphi_{n}(u)w_1(u)du\right]w_2(v)dv = 0 \\
  \implies \int_{a}^{b}f(u,v)\varphi_{n}(u)w_1(u)du = 0 \mbox{ a.e. $v$ } \\
  \implies f(u,v) = 0 \mbox{ a.e. $u$, $v$ }.
$$
There are details to be worried about concerning a.e., but they generally work out just fine because the eigenfunctions are smooth. The functions $\{ \varphi_n(u)\psi_{n}(v) \}_{n,m}$ then form a complete orthonormal basis of eigenfunctions in the $2$-d. case, which means you have a $2$-d expansion
$$
        f(u,v) = \sum_{n}\sum_{n}(f,\varphi_n\psi_m)_{w_1w_2}\varphi_n(u)\psi_m(v).
$$
You can extend these arguments to higher dimensions; there are no real gotchas. And, I think you can probably see how to bootstrap to higher dimensions. Once you know that everything can be expanded using the separated eigenfunctions, then there can't be any more. It's an interesting argument that's not too difficult, and I'll leave to you: If you have an orthonormal basis of eigenfunctions of a selfadjoint operator, then there cannot be any others. (Hint: expand the new eigenfunction in terms of the others.)
The infinite intervals are trickier because the expansions may include discrete spectral elements and continuous ones. For example, on an interval $[a,\infty)$, the expansions may be
$$
              f = \int_{\sigma}\left(\int_{a}^{\infty}f(t)\varphi_{\lambda}(t)w(t)dt\right)\varphi_{\lambda}(x)d\rho(\lambda),
$$
where $\rho$ is a spectral density measure on the spectrum $\sigma$ of the operator, and $\varphi_{\lambda}$ are eigenfunctions with real eigenvalues satisfying a condition at $a$. The eigenfunctions may not be square integrable,
which means they're not really in the space, even though integrals in $\lambda$ over a small range of eigenvalues give wave packets that are square integrable, and are nearly eigenfuntions. The above may have discrete components and or continuous components, but the idea that you can expand somehow in the eigenfunctions is still the same. The are Parseval identities for such expansions, and properties that make these continuous and/or discrete expansions look like the Dirac formalism. (Keep in mind that Dirac would have have known about such expansions from Weyl, Hilbert and others.)
The cases grow in complexity with singularities, but the basic conclusion is this: you'll probably never have to worry about things not working the way you expect in regard to completeness. :)
