Let our domain $\Omega = (0,2\pi)\times(0,2\pi)$ be a square with sides parallel to the axes.
Consider the Dirichlet eigenvalue problem $\Delta u+\lambda u=0$ with Dirichlet boundary condition $u=0$ on $\partial\Omega$. It is known that separation of variables yields the solutions, for example in these lecture notes.
My question: Why are all solutions given this way? Why doesn't there exist a solution that can't be obtained by separation of variables?
(Of course, a linear combination of solutions obtained by separation of variables might not have separated variables, but I'm asking why there isn't something other than those basic solutions and their linear combinations.)
For a partial differential operator $L$ with a domain $D_L\,$, this is always the case whenever $L$ admits separation of variables on $D_L$. The variables are to be separated successively one-by-one. Your example of $L=\Delta$ on $D_L=H^1_0(\Omega)\cap H^2(\Omega)$, with $\Omega\subset\mathbb{R}^2$ being a square, is just the simplest instance of such successive separation of variables. It is absolutely unimportant what variable shall one start with. Let it be $x$. In this case, consider the Sturm-Liouville problem $$ \begin{cases} X''=\mu X,\;x\in (0,2\pi),\\ X(0)=X(2\pi)=0 \end{cases} \quad \Rightarrow\quad \begin{cases} X_j(x)=\sin{\frac{jx}{2}}\,,\\ \mu_j=-\frac{j^2}{4}\,,\;j\geqslant 1. \end{cases} $$ Due to the widely known fact that $\{X_j\}_{j=1}^{\infty}$ is an orthogonal basis in $L^2(0,2\pi)$, every eigenfunction $u\in H^1_0(\Omega)\cap H^2(\Omega)$, i.e., any nontrivial solution $Lu+\lambda u=0,\;u\in D_L\,$, can be developed into Fourier series $$ u(x,y)=\sum\limits_{j=1}^{\infty}Y_j(y)X_j(x), \quad Y_j(y)=\frac{(u,X_j)}{\|X_j\|^2} =\frac{1}{\pi}\int\limits_0^{2\pi}u(x,y)\sin{\frac{jx}{2}}dx,\;\,j\geqslant 1. $$ Applying the habitual multiplying and integrating by parts, one readily gets BVPs for the Fourier coefficints $Y_j\,$, namely, $$ \begin{cases} Y_j''+\Bigl(\lambda-\frac{j^2}{4}\Bigr)Y_j=0,\;y\in (0,2\pi),\\ Y_j(0)=Y_j(2\pi)=0,\quad j\geqslant 1. \end{cases} \tag{1} $$ It is clear that nontrivial solutions $Y_j$ of BVP $(1)$ exist if only $$ \lambda-\frac{j^2}{4}=\frac{k^2}{4}\,,\quad j\geqslant 1,\;k\geqslant 1. $$ Hence every nontrivial solution $Lu+\lambda u=0,\;u\in D_L\,$, is of the form $$ u_{jk}(x,y)=\sin{\frac{jx}{2}}\sin{\frac{ky}{2}}, \quad \lambda_{jk}=\frac{j^2}{4}+\frac{k^2}{4}\,,\;\;\,j,k\geqslant 1, $$ except for certain eigenvalues $\lambda$ which admit multiple representation $$ -4\lambda=j_1^2+k_1^2=\dots=j_N^2+k_N^2,\quad N=N(\lambda)\geqslant 1,$$ where $N(\lambda)$ denotes a geometric multiplicity of $\lambda$. Here is the list of the first 24 multiple eigenvalues, i.e., with geometric multiplicities $N(\lambda)\geqslant 2\,$: $$ \begin{align} -4\lambda=50 = 1^2 + 7^2 = 5^2 + 5^2,\\ -4\lambda=65 = 1^2 + 8^2 = 4^2 + 7^2,\\ -4\lambda=85 = 2^2 + 9^2 = 6^2 + 7^2,\\ -4\lambda=125 = 2^2 + 11^2 = 5^2 + 10^2,\\ -4\lambda=130 = 3^2 + 11^2 = 7^2 + 9^2,\\ -4\lambda=145 = 1^2 + 12^2 = 8^2 + 9^2,\\ -4\lambda=170 = 1^2 + 13^2 = 7^2 + 11^2,\\ -4\lambda=185 = 4^2 + 13^2 = 8^2 + 11^2,\\ -4\lambda=200 = 2^2 + 14^2 = 10^2 + 10^2,\\ -4\lambda=205 = 3^2 + 14^2 = 6^2 + 13^2,\\ -4\lambda=221 = 5^2 + 14^2 = 10^2 + 11^2,\\ -4\lambda=250 = 5^2 + 15^2 = 9^2 + 13^2,\\ -4\lambda=260 = 2^2 + 16^2 = 8^2 + 14^2,\\ -4\lambda=265 = 3^2 + 16^2 = 11^2 + 12^2,\\ -4\lambda=290 = 1^2 + 17^2 = 11^2 + 13^2,\\ -4\lambda=305 = 4^2 + 17^2 = 7^2 + 16^2,\\ -4\lambda=325 = 1^2 + 18^2 = 6^2 + 17^2 = 10^2 + 15^2,\\ -4\lambda=338 = 7^2 + 17^2 = 13^2 + 13^2,\\ -4\lambda=340 = 4^2 + 18^2 = 12^2 + 14^2,\\ -4\lambda=365 = 2^2 + 19^2 = 13^2 + 14^2,\\ -4\lambda=370 = 3^2 + 19^2 = 9^2 + 17^2,\\ -4\lambda=377 = 4^2 + 19^2 = 11^2 + 16^2,\\ -4\lambda=410 = 7^2 + 19^2 = 11^2 + 17^2,\\ -4\lambda=425 = 5^2 + 20^2 = 8^2 + 19^2 = 13^2 + 16^2. \end{align} $$ Since all eigenfunctions $u_{jk}$ are pairwise orthogonal, for every particular eigenvalue $\lambda$, a Fourier series for solution corresponding to this $\lambda$ degenerates into a finite sum $$ u(x,y)=\sum_{r=1}^{N(\lambda)}c_r(j,k)u_{jk}(x,y). $$